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The final chapter of this report pertains to a FORTRAN computer 
program which implements and numerically integrates the complete set of 
equations* The salient features of the program, its subroutines, and 
the input and output data are described. An annotated flowchart along 
with a full listing of the code is provided. 

It is noteworthy that while the idealization and methodology ap- 
plied in this report are essentially those of Likins, ^ ^ the equations 
formulated herein are unique from those developed in Reference 5, and 
indeed the distinction is fundamental- It was the express desire to 
avoid the kinematic restrictions required there to effect a coordinate 
transformation on the elastic deflections which motivated this approach. 

From an applications standpoint, the basic discretization of the 
vehicle of interest is performed in the manner of lumped mass structural 
dynamics modeling. The required stiffness matrix which quantifies the 
internal elastic restoring forces can in general be obtained from pre- 
processed linear structural finite element analysis programs (e.g., 
NASTRAN). Because of the mass particle idealization of the elastic do- 
main, only translational displacements are defined at those points, 
hence any finite element model used to provide stiffness matrix informa- 
tion must be purged of any rotational degrees of freedom that may exist. 
This requirement is easily satisfied through the application of the 
static condensation procedure. Thus the analyst is afforded these fa- 
miliar and versatile structural modeling techniques augmented by the 
arbitrary motion capability. . 

The motion equations formulated here are complete and unabridged 
for a single unconstrained flexible vehicle. However, they could, in a 
straightforward fashion, be coupled to the dynamics equations for other 
independent bodies to form an articulated system. This can be done 
through the identification and elimination of interbody constraint forces/ 
torques and redundant kinematic variables. Indeed, it is for just such 
an application that these equations are intended. Specifically they are 
to represent a generic flexible payload to be terminally attached to the 
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Space Shuttle Orbiter remote manipulator system, which is an articulated 
chain of rigid and flexible bodies. For this case the model's rigid-body 
is taken to be the payload grapple fixture with all outboard structure 
represented by the particle assemblage - 
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CHAPTER 2 


PRELIMINARIES 


2 . 1 Vehicle Idealization 

The system being analyzed (see Figure 1) consists of a single 
rigid body and an attached flexible appendage. The appendage is ideal- 
ized as a system of particles connected by massless elastic structure. 

There is no articulation between the appendage and rigid base, i.e., 
the appendage is "cantilevered" to the rigid body. At an arbitrary 
point, 0 , of the rigid body we locate the origin of the body fixed frame 
which rotates as the body rotates in inertial space. The vector R serves 
to determine the position of Og* relative to the inertially fixed point 0. 
The particle masses m^ (i = l,2,...,n) are located via the position vectors 
r. relative to Og in the unde formed state. The elastic displacement of 
m^ is measured in the body frame. 

Many space vehicles or parts of spacecraft can be approximated in 
this manner. A specific example is the Shuttle Remote Manipulation Sys- 
tem in which the "appendage" corresponds to a flexible payload and the 
"rigid base" to the grapple fixture (this component being attached to 
the orbiter through the links of the manipulator arm) . 

2.2 External and Internal Forces 

With the ultimate goal in mind of applying the present analysis 
to more complicated situations , we wish to accommodate all forces and 
torques which will arise when the system in Figure 1 is attached to other 
spacecraft components. Hence, at the point Og, let there be a force-torque 
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pair: f (t) , T (t) . In the domain of the elastic appendage we have an 

external force f (t) acting upon the point mass m^. 

As an example, in the case of the Shuttle Remote Manipulator Arm, 
f and T would represent the force and torque exerted by the end-effector 
on the grapple fixtxire. 

We assume that we are given a stiffness matrix reflecting the 
mutual elastic forces between the mass points of the appendage. Assemble 
the elastic displacements as 

,111222 nnn,T 

3 = (q^ qy qy 

(x,y, 2 ) refer to Cartesian components along the axes of the body fixed 
frame . 

If [K] is the stiffness matrix, then -[K]q is the vector of elastic 
forces exerted on the point masses. Assxime that [K] is partitioned such 
that the vector of elastic forces are ordered exactly as the elements in 
Note that in generating [K] the appendage is cantilevered to the 
rigid base. 



Figure 1. Idealized vehicle. 
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CHAPTER 3 


EQUATIONS OF MOTION 


3.1 Particle Translational Equations 


. th 


V t the inertial velocity of the i particle, is given by 


~^i d ”^i -»“i , 

V = — (R + r+ q) 

at 


ru> 


Let I y ) and ^ represent the (absolute) velocity and angular velocity 
of the body frame resolved in body axes, we then have 


i #1 , 1 I, *1 


w> 


Differentiating this expression we arrive at the particle acceleration 


u 


^ i (*1 /i. ^ ^ ^ 

— V =lvj~(r +£)xw + 3 
\w/ 


+ (jJ X 


V j + 2a 

VWi 


+ ^ X X + a^) ] 
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Expanding the cross product in the last term and using the matrix-vector 


form for the cross product ^ ^ h - [al^b where [ £ ] " = a 0 


° -^3 ^2 


-a a 

L 2 1 


the particle acceleration can be written as 


d i 
dt - 


V - ([r^]~ + + [w]~ V + 2 [w]~£^ 


,Ti. 1!I|2 i ,Ti. I |2i 

“ 1 llilM £ + (w £ “ 1 1 £ 


(3-1) 


If f^ is the external force on the i particle and ^ the elastic force 
exerted on the r particle by the rest of the assemblage (both resolved 
along body axes ) , the translational equation is 


d 1 

f + f = m. -r— V 

— -e 1 dt — 


(i “• lf2/**«fn) 

Partitioning the stiffness matrix into (3x3) arrays 


[K] = / [K 





“^21> 


■ • 

“'nl' 

[k'jI.. 


n 

E 

[K, .1 


j=i 



Employing Eq. (3-1) for the particle acceleration, the translation equa- 
tions for the appendage particles may be assembled as 
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m 


in^{r + 3 ) 

~ 2 -2 

+ £ ) 


Lm^ (r + 3 ) J 


0) + 


m 

0 

0 


0 0 

2 ^ 

m 0^ 

0 m'* 


0 

0 


+ [K]2 = 


+ u 


(3-2) 


where we have introduced the symbol m 


0 m. 0 

1 

0 0 m, 

1 


(i — l/2/.*>fn) 


and the nonlinear kinematic term u is given by 


Hv 


— 


- • 1 




m^(u 3 



/u \ 




/ \ 

- *2 


"^2— 

w 



• 

\ / 

. 


• 


• 





• 


Z 


- *n 


m ^ 


m ^ q 


_ n M 


^ n — 



/I lx “ 

m^^ • (r + £ 

2 2 

m ^ * (£ + 3 )w. 


/ n n. 
m ^ • (r + 2 ) 


+ 0) 


" . 1 ir 

m^ (r + £ ) 


2 2 
m (r + 3 ) 


, n n. 
m (r + q ) 
- n — 


(3-3) 


Equation. (3-2) constitutes a set of 3n scalar differential equations. 


9 


3.2 


Vehicle Translational Equations 


For the composite system (rigid body and appendage) the sum of 
the external forces equals the total mass times the acceleration of the 
mass center. 

If is the mass of the rigid body and _s is the vector from 
to the mass center of the rigid body (expressed in the body frame) 

n 

m£ = (£. + 3 ) + (3-4) 

i=l 


n 

m = m^ + ^ m^ is the total mass. 
i=l 

c(t) is the vector position of the instantaneous mass center relative 
to 0 . 

d 

The acceleration of. the mass center is: — 2 (R + c) 

dt 



dt 



(in body frame) 


m 



n .. . ^ n ^ , 

m.q^-m£x^+ 2 ^x ^ (^x me) 

i=l ^ i=l 


2 

Expressing this last term as: (£ • m£)£ - | l^l [ m£ the vehicle trans- 

lational equation assumes the form 


Im 0 o1 


l-\ 

1 


• 

o 

B 

o 


V 




B 

o 

S. 


\ w / 

— j 



[mej^o) + [m 


m 


m ]£ 


E £ 


-f 


i=0 


(3-5) 
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The nonlinear term u 


given by 





n 

2b} ^ m.q 

- i=i ^ 


( oi • me ) oj + 



(3-6) 


3. 3 Vehicle Rotational Equations 

For the composite system (rigid body and appendage) the s\im of 
the external torques taken about the mass center equals the time rate of 
change of the angular momentum taken about the mass center. 

Let [I, ] be the inertia matrix of the rigid body with respect to 
b 

a coordinate system located at the mass center of the rigid body and 
parallel to the body fixed axes system at 0^. 





Xe - xx"^] dm 


^ is the position vector of a mass element dm in the rigid body relative 
to the rigid body mass center and the integration is performed over the 
region occupied by the rigid base, [E] denotes the unit matrix. 

The system angular momentum can be split into two parts: 


center 


4 


— angular momentum of rigid base 


n 



Let ^ 

to m. 

1 


— angular momentum of appendage particles 

and ^ denote the position vectors from the system mass 
cind a generic mass element in the rigid body respectively. 


4 ■ 


H. = X m. ^ 

-i — 1 dt — 
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Figure 2. Vector geometry. 


From Figure 2, ^ A ~ A “ £• 

Inserting this expression for £ into the integral definition of 
^ and recalling the definition of and the fact that ffh dm = 0_ 
we arrive at the following expression for ^ 


4 


[i^]^ + 




c) 



(s - c) 


Thus 




[I^]w + £ ^ [w X £ + W X X £)] 


m^ (£ - £) 
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Turning to the angular momentum of the i particle 


at 




?L “ *^3^) X ^ + 2^ X + [^ • (r^ + £^)]^ 


j I 1 2 i i d 

kl I (£ + a ) — 


Combining the above expressions for and (with the substitu- 

d2 ^2 

tion for — - ^ ) the terms involving — ^ £ conveniently cancel leaving 
dt"^ dt 

the following result for the time derivative of the angular moment lam 


II 

'it- ^ - m^(£ - c)~^)w + w X ^ 


- c) X [u X (w X s^)] - £ *11^ tl!*’] ' + 3^)'“ 


+ 2 £ m. ^ [w. • ^ 

i=l ^ i=l ^ 


1 ^ 11 ^ 2 ^ + 3 ^) 

i=l ^ 


(3-7) 
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The system rotational motion equation is 


^ H = 1° - £ ^ f ° + 2 

at 

= 1° - £ ^ X] 3^) 

j=0 i=l 

If we were to insert Eq. (3-7) into this last equation we would have a 
valid equation for system rotation. Note, however, that this equation 
would not depend upon (u,v,w) explicitly- Since the translational equa- 
tions (3-5) depend upon ^ explicitly, and we wish to have a final set of 
equations with a symmetric coefficient matrix of the generalized accel- 
erations, we force the coupling between the rotational equations and the 

• • • 

acceleration vector (u v w) by the following device. 

n 

Take Eq. (3-5) and (3-6) and solve for ^ The system rota- 

j=0 

tional motion equation then becomes 


u\ 


f H 
dt - 


+ 


i=li 


+ [c] [me] 0 ) 


/ 


- £ X ^ m . - [ m £] *"^ vj - 2[£]‘'£_2 m .£ 

i=l* ^ \w/ i=l ^ 


(tj • me) c X 0 ) 


(3-8) 


Combining Eq- (3-7) and (3-8), we arrive at the final desired form for 
the vehicle rotational equation 
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[n>£] ~ 



[I(t)] 


— [j"! 


1 II 2 2 ' 

+ £ ) ' I m 2 (r +£)']. 


I , n 
I ni (r 
I n — 



o 

= T + 


Z 

i=l 


/I 

(r + 3 ) 




(3-9) 


Here 


n ... 

[I(t)] = [Ij^3 - mj^(s_ - £)~i. - Z m^U^]~(r^ + £^)~ - [£] [m£] ~ 

i=l 


which can be simplified to 


[I(t)] 


r 1 r ~ 1 2 V' 2 

^ ~ "^3 - ~ Z "‘i - ^ ^ 

i=l 


(3-10) 


In this form we recognize [I(t)] as the inertia matrix of the vehicle 
about point 0 . 

g 

The nonlinear rotation term u is given by 


u 

— r 



n 

2[£]"'[^]"' 2 ~ * m£) [£] 

i-1 ^ 


n 

- “ ^ (£ “ £) ^ [£ X (£ X ^) ] - 2 f 

i=l 


11 ... ^ 11 
- ^ [£ • (£^ + £^)]m, ^ X £ + I jul I ^ m. X (£^ + £^) 

i=l ^ i=l ^ 


(3-11) 
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CHAPTER 4 


EXPANSION AND PARTITIONING OF TERMS 


4. 1 Deformation Dependent Coefficients 

In the equations developed thus far, specifically Eq. (3-2), 

(3-5) , and (3-9) , we have isolated the accelerations on the left hand sides 
of the respective equations. The acceleration coefficients are time 
dependent through the elastic deformations. It is quite desirable from 
an applications viewpoint to rank the constituents in these coefficients 
in accordance with their relative magnitude. Thus, in a computer simula- 
tion, one can choose to omit certain terms and speed up execution with a 
minimal impact on computed results. 

We will rank terms amongst three categories : 

(1) Terms independent of 

(2) Terms first order in 

(3) Terms second order in £. 

The majority of coefficients are directly identifiable in this hierarchy. 

We have two coefficients which require additional attention: m£ and 

[I(t)]. 

n . n . 

From Eq. (3-4), me = itL £ + ^ first 

^ i=l ^ i=l ^ 

two terms are of category 1 and the sum will be denoted by The 

time -dependent term will be denoted by mc^(t). 

The matrix [I(t)] is given by Eq. (3-10) and can be written as 
[I(t)] = [1^1 + tl^Ct)] + [l3(t)l 
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The three matrices and [I^(t)] are of category 1, 2 , and 

3 respectively and are given by 


-i, 2 


[111 = [Ij^] - 


[I (t)] = m {[r^][ 2 ^] + [i^][£^]) 

i=l 


~i ,2 


[I (t)] = m ] 

^ i=l 


4 . 2 Nonlinear Kinematic Terms 


In this section, we concentrate upon the three nonlinear terms: 

u , u » and u appearing on the right hand sides of the motion equations. 
— t ~v ^ 

Following a procedure similar to that of the previous section, the non- 
linear terms are partitioned amongst three categories: 


( 1 ) 

(2) 

(3) 

Accordingly, 


Nonlinear terms independent of 3 , £ 
Nonlinear terms first order in £, £ 
Nonlinear terms second order in £, £ 
from Eq. (3-6) , u^ = u^ ^ + u^ ^ with 




■** 



(4-1) 



n 

= m 

i=l 


• Itl£^)^ + 



(4-2) 
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In a similar manner from Eq. (3-3) 



■2 V m. ([c]' + Li"-]') = 

i=l ^ ■ 


-2^ m. (r^ + 
i=l ^ 


n . ^ . n 

= -2^ m. [r^] [w] - 2 y^m [g^] " M 

i=l ^ i=l 
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The third term in can be expressed as 

(u * m£)[c]~(;). = (W • + (u • m^) + (w * mc^) 

+ • in£^) [£^] 

For the last two terms in the following expansions are useful 


• (r^ + ^ 


(^ • (£^ ~ ^ — 




"*■ " £l^ ^ £ 




X (r^ + £^) = X ^ + r"- X + 3^ X Cq + s'" ^ £l 


Collecting terms in u independent of deformation 


( 1 ) 

u 

— r 


j 


[m^]*' [£]“* I V j - (£ • mc^) [c^]'“£ - 


(s - Cj^) X [£ X (£ X £) ] “ ^ '! ^ ~ ^ ^ — 

— “0 • ^ 


i=l 


+ 1 1“M^ £ ^ 

i=l 




The expression for can be further simplified by use of the follow- 

ing identities 
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£ "'i ^ 

i=l 






5^ m . (^ • ) (r^ " £q ) ^ iJi = “ tiill ~ £.^ £,^ 

i=l i=l 

+ [(£ • (nij^ - mc^) ] Cq X w 


-nij^(£ - ^) X [^ X (u X £) ] = - ( w • £)n>j^ ^’'•fL + 

- 1 kl 1^ \ ^ E 


Incorporating these results we arrive at the final expression 


( 1 ) 

u 

—r 


u\ n i iT 

- [mc^] [w] "* V + [£] ~ 2 m. r^ ^ “ il ^ — 

-0 \ i i^l 1 


i wy 


( s_ • ^)_s X ^ 


( 4 - 5 ) 


Collecting first order deformation dependent terms in u 


( 2 ). 


' u\ n 


-[mc^]~[w]~ |v - 2^ m^[£ ]~[w ]~3 - (u • mcQ)[c^]' 

\w/ i=l ^ 


n . . ^ ^ 

[ (^ • (£^ ” ^ ^ ^ ^ + (^ • £. ) (3 ” — 1 ^ ^ 


i=l 


+ 1 kll^£ ^ El + E ^ £^'> 

i=l 
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The expression for can be further simplified by use of the follow- 

ing identities 


^ i 

• _*i ■*■ 


"E X W + (<»).• £^)n»^(a^ - ^ = 


i=l 


[w]~E (£V^ + + (£ • ”‘£i^£o "" - 

i=l ^ 




2 

X [^ X X £) ] = • £)£^ ^ 1 S.^. ^ — 


( 2 ) 


/u\ n ^ 

-[me ]'[u]~ V - 2^ m^[r ] ' [m] ~£i 
\w/ i=l 


_ , i iT i iT. 

+ [w] 2-( (£ a + a £ )£ 

i=l 


( 4 - 6 ) 


Collecting second order deformation dependent terms in and simplifying 


( 3 ) 


^ *i r i iT 

= -2V m.[q. ]'[(»)] £+[£]>, m £ a (o 

i^l ^ itl ^ 


( 4 - 7 ) 


u 

— r 
Eq. 


where the terms on the right hand side are given by 

— r — r — r 

( 4 - 5 ) - ( 4 - 7 ). 
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CHAPTER 5 


MODAL COORDINATE TRANSFORMATION 


When the number of particles in the appendage idealization becomes 
large, high frequencies obtain which make numerical integration difficult. 
We will describe a truncated coordinate transformation to circumvent this 
difficulty. Note that the treatment to follow is somewhat heuristic and 
hence requires good engineering judgement and caution in its implementation 

Since the high frequencies arise from the appendage vibration, it 

is natural to start with its governing equation (3“2, 3“*3). Consider the 

• • • 

case where no external forces act, ^ and (u,v,w) = The ’‘constrained 

appendage equation then assumes the familiar form 

[M]£ + [K]£ = 0 where [M] = diag (m^ ,m^ , . . . ,m ) 

The natural frequencies, , and corresponding mode shapes, v^, are 
determined from 

([K] - w^[M])V = 0 (5-1) 

For the vehicle we are treating here, the appendage is rigidly 
attached to the base body hence no rigid body modes are present in the 
above eigenvalue problem- Equivalently, [K] is positive definite- Since 

[K] and [M] are symmetric and positive definite there exists 3n independ- 

i 2 

ent eigenvectors v corresponding to positive eigenvalues o)^, even if 

there are multiple eigenvalues. 
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We assume that the eigenvectors are normalized such that (v^/[M]v^) = 1. 

It follows that (v^,[K]v^) = We can always create a mutually ortho- 

— — —1 

gonal set such that 

(v^,[M]v^) = 0 = (v^,[K]v^) (i j) 

In actual computation we can deal with a simpler eigenvalue problem 
than that presented by Eq. (5-1). Specifically we will transform Eq. (5-1) 
to an ordinary symmetric eigenvalue problem. Introduce the change of 
variables; W = Since . [M] is diagonal, [M] is a diagonal 

matrix whose elements are the square roots of the corresponding elements 
in [M] . The eigenvalue problem transforms into 

or [^]W = o)\ (5-2) 

- 1/2 - 1/2 

where is the symmetric matrix: [^] = [M] [K] [M] 

It is easily verified that if the eigenvectors W^ (i = l,2,...,3n) 
of Eq. (5-2) are orthogonal (which can always be done) then the corre- 
sponding eigenvectors of (5-1) satisfy all orthogonality and normality 
conditions specified above. 

2 2 2 

Order the eigenvalues such that < o)_ < ... < oj and let [^>] 

^ ^ ^12 t 

be the (3n x t) matrix whose columns are the eigenvectors v ,v , . . . , v 
(t ^ 3n) . We now make the transformation 

£ = l^ljl (5-3) 

This is not a coordinate transformation in the strict sense, since [^] 
does not have an inverse when t < 3n. The appendage deformation is now 
characterized by t ' "modal coordinates" instead of the original 3n defor- 
mation coordinates. We formally make the substitution, Eq. (5-3) , into 
the full set of motion equations. 
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Substituting into the appendage defcyrmation equation (3-2) , pre- 
multiplying by and recalling the orthogonality and normality con- 
ditions we arrive at 


: V 


/ 1 1 .^ 

(r + £ ) 

2 2 ^ 

+ a ) 


/ n n_ 
m (r + q ) 

- n — 


(jj t n + 


T T 

[$] . + [^] U, 


(5-4) 


The vehicle translational equations (3-5) become 


m 0 0 / u\ 

0 m 0 I V - [m£] + [m^ m ... m^] [^]ri 

00m \w/ 


i=0 


The vehicle rotational equations (3-9) become 


me ] " V + [ 


• r 1 1 • .2 2, ^ > 1 , n n, , . *• 

;i(t)]£ + 1^2 ~ ^ ^ “ 'i "'n — ^ ^ J ^ ” 


^0 + E 

- i=l 


The assembled equations of motion in matrix form are presented in Figure 3. 
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CHAPTER 6 


SYSTEM KINETIC ENERGY 


The kinetic energy of the vehicle is the sum of the translational 
and rotational kinetic energy of the rigid body and the kinetic energy 
of the particles comprising the appendage 


1 2 1 , , 1 ^ 2 
= I ^ \ ^ I ^ I “i '^i 


V = (v| + a)Xsis the velocity of the mass center of the base 

W/ ' “ 


i 

V = V 


W/ 


* • * th 

+ 0 ) X (r^ + 3 ^) +3 is the velocity of i particle 


Forming the inner products recalling Eq. (3-10) for 

[I(t)l the kinetic energy can be written as 


, n . ™ 

1 o O 2 1 T 1 T ' ^ • i X *1. 

= *rin(u + v + — ^[I(t)]^+ — 


i=l 


1 1 T I i '1 

- — (uvw) [mc]"'^ + T ~ "9 

2 W/ ^ i=i ^ 


+ — 


n 

(“1 

rr. ^ 

T n 


1 1 T V' / ^ i 1.^*1 

+ 7 U) (r + 3)3 

-IS 

1 — 1 
n 

•H 


i=l 

i=l 
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• i • 

We now rewrite those terras in T which depend upon £ in terras of r[ 


E .-Ini**; *rn • »TT • *T • 

{^ } m.q = 3 [M]q = [M] [$] jn = n[E]n 

i=i ^ 


[E] is the (t ^ t) identity matrix 


n 

(uvw) ^ 
i=l 


m. q 


(uvw) 


r I'l 2 

m [ m 


(uvw) [ra' 


2 

m 


•1 



. [$] n 


m^(r^ +£^)' ••• + a")~] [4>]n 

i=l 

1 T 

The kinetic energy cam be written as the quadratic form T = — U [A]U 
where [A] is the coefficient matrix (symmetric) of the generalized 
accelerations appearing in the equations of motion (see Figure 3) and U 
is the vector of non-holonomic velocities 

/ I T * *T \ T 
U = \uvwj^ \x]_ j 


Since [I ] is positive definite, an inspection of the initial 
expression for T reveals that T ^ 0 for all JU. If T = 0 then v^ “ ^ ” 

v^ = 0_ (i = 1,2,:.’. ,h') . But == £ = iiL iu^plies Cuvw) = _0 and = £ = 


(uvw) = ^ implies q = Hence [$]n_ = £. Since the columns of [<I>] are 
linearly independent we > must have = £ also. In other words, T = 0 if 
and only if £ = £• This argument proves that [A] is positive definite 
and consequently nonsingular (see Chapter 8 where we require [A] ) . 
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Note that if we replace the rigid body by a particle then ^ = 0^ 
and [Ij^] = [0] . We still have T ^ 0 but if T = 0 we can only argue that 
(uvw) =0. We can have T = 0 for nonzero u and _n as long as w x 
+ = 0 (i = 1,2 , . . . ,n) . Thus for this later case [A] is positive semi- 

definite. In particlar [A] will be singular. The situation here can be 
understood by simply enumerating the degrees of freedom involved. Orig- 
inally we had a system consisting of a rigid body and n particles: (6 + 3n) 

degrees of freedom. The number of dynamic equations was also (6 + 3n) . 

When degenerating the rigid body to a particle we have a system of (n + 1) 
particles: (3n + 3) degrees of freedom. However, when we retain the same 

equations of motion as in the original case (6 + 3n) there will clearly be 
a redundancy present. Indeed, this explains why [A] is singular for the 
degenerate case. Consequently, we cannot use the equations developed here 
for a system composed solely of particles; at least not without modification. 
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CHAPTER 7 


KINEMATICAL RELATIONSHIPS 


Let the transformation from the inertial frame to the 

body frame {x ,y ,z } be arrived at by a sequence of three Euler angles 
depicted below. 




12 


34 


R 



0 


0 

\ 

1 

/ cos 6^ 0 

sin 

0 

cos 

1 

1 — 1 

CD 

•sin 

1 — 1 

CD 

= 

0 1 

0 


sin 


cos 



pin ^2 ^ 

cos Q^i 

^cos 

S 

-sin 



i 



sin 


cos 


0 

[R^^] is the transformation 

0 


0 


li 

f matrix from frame ’ j * 

to 


frame ' i * 


14 12 23 34 

Concatenating transformations, [R ] = [R ] [R ] [R ] 
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[R 


14, 


i cosS^cosB^ 

cos0^sin9^ + sinS^sinO^cosO^ 
^sin6^sin6^ - cosS^sinS^cosO^ 


-cosO^sinS ^ 

cos6^cos9^ - sinO^sinS^sin©^ 
sin0^cos0^ + ccs0^sin02sin0^ 


sin 62 

-sin0^cos02 

cosO^cos©^ 


(7-1) 


We next derive the relationship between the body frame angular 

• • # 

velocity to (expressed in body coordinates) and the Euler rates ^2' ^3* 

Let {i ,j ,k } be the set of unit vectors along the axes of frame 'p' 

. -P -P 
(p = 1,2, 3, 4) 


= ®lH 


• 

+ 0-,j 


2-'2 ®3^^3 


To express uj in the body frame f we must use the representation of 
the unit vectors in frame 4. With the aid of the transformations listed 
above we arrive at 


/• • \ 

/ 0^cos02COs 9^ + S^sin©^' 

©^cos©^ - e^cos©2sin©^ 

©3 + e^sin©2 


(in body frame) 


This system can be inverted to yield 


n \ 

1 

2 

\“3/ 


cos 9. 
cos 9, 


sin 0. 


-cos © 3 tan ©^ 


-sin 0. 
cos 0, 


cos 9, 


sin © 3 tan ©^ 


(cos ©3 ^ 0) 


(7-2) 
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CHAPTER 8 


EQUATIONS OF MOTION FIRST ORDER FORM 


The assembled motion equations (Figure 3) can be written as 


where 


[A(t)] 


u 


0 




V 


0 





w 

= F + U - 




1 1 



2 

« 





2 2 

n 





2 


0) ri^ 

1 1 

_t tj 


( 8 - 1 ) 


F = 



(8-2) 
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and 


U = 


u 


T 

u 


(8-3) 


Let R f R f components of the inertial position vector o 

(origin of body frame) resolved along inertial axes and [V] denote th 
matrix in Eg. (7-2) . The )cinematic relationships can now be written 


lk\ 


R 


= [R^^l I V 





= [r]0) 


Define [Q 1 


= diag (ui^coi - . . • The state vector Y is defined to be 

12 r. 


T T *T,T 

Y = (R R R ^ ■!! ^ 

_ X y z 1 2 J— 


(8-4) 


The equations of motion written in 


first order form are 


dt 


Y 





[T]i^ 



(8-5) 


This system of (2t + 12) first order equations can be integrated numer 
ically with appropriate initial conditions. 
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CHAPTER 9 


DIGITAL SIMULATION 


This chapter is concerned with the FORTRAN computer program which 
implements and numerically integrates the complete set of first order 
ordinary differential equations presented in Chapter 8 , Eq. (8-5). A 
description of the main program, its subroutines, and the input data is 
given. An annotated flowchart of the program is given in Figure 4 and 
a complete listing of the program and its sxibroutines is provided in 
Appendix A. An example of the input data for a sample vehicle is pro- 
vided in Appendix B. The code is liberally commented throughout and in 
most instances the FORTRAN variable names are mnemonically similar to the 
corresponding analytical quantities. Virtually all computations involving 

real nijmber quantities are performed in (IBM) do\i>le precision. External 

( 8 ) 

subroutinBs front th© doulDle prscision IMSL libiraxy are used, to perform 
certain standard computations. IMSL subroutine ”EIGRS is used for 
03 _g 0 nvalue /eigenvector extraction and subroutine LEQTIP is used to 
solve simult 3 ^neous linear equations. In addition, IMSL subroutine USPLT 
is used to generate time history graphs of selected elements of the vector 


T T 

{r R R 6, q ••• vivw OJ oj w ^ ... 2 } (9-1) 

xyzl23-* ^ izj, 

via the line printer. 

Throughout the program deformation dependent terms are arranged and 
computed hierarchically as quantities involving structural deflections to 
the first and second degree. Similarly the nonlinear kinematic terms 
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START y 

^ziz: 


SET INITIAL TIME 


SUBROUTINE INITL 


SUBROUTINE GNMASS 



SET APPLIED FORCES 
AND TORQUE 


t 

SUBROUTINE EXTF 



SUBROUTINE SOLVE 





UPDATE TIME 


READS INPUT DATA 

COMPUTES EIGENVALUES & E'GENV|CTORS 
COMPUTES TIME INVARIANT PARAMETERS 

COMPUTES & ASSEMBLES DEFORMATION 
dependent GENERALIZED MASS MATRIX 


COMPUTES & ASSEMBLES VECTOR OF APPLIED FORCES 
assembles & INTEGRATES STATE EQUATIONS 



PRINTS TIME. APPLIED FORCES, 
ySUBROUTINE PRIlTry KINEMATIC VARIABLES. 

Z_ ^ ' AND STRUCTURAL 

deformations 


/subroutine PLOTy / PLOTS SELECTED VARIABLES VIA LINE PRINTER 


C 


STOP 


3 


Figure 4 . Program flowchart. 
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are organized into the three categories of Section 4.2 with the contribu- 
tions of each group of terms being computed independently. This parti- 
tioned structure of the computations provides the capability to assess 
the influence of these higher order terms on the final solution and upon 
such analysis bypass those deemed negligible. 

Main Program 

The main program is simply an executive module which calls the 
appropriate siabroutines in the proper order. The reader will note ^that 
if external forces are required, these must be explicitly coded either 
in the main program or as individual sxibroutines . If the external forces 
are time dependent, it is essential that they be recomputed prior to each 
call to subroutine ”EXTF" (see comments in main program) . For the sys- 
tem in Fig\are 1 the external excitation is accommodated via the three 
* 

arrays: F(?, TAU(2f, FP. 

F(2f — s\am of external forces on rigid body 

TAU? — sum of external moments on rigid body taken about body 
frame origin. 

F0 and TAUflf are three-dimensional vectors whose elements refer to com- 
ponents along body frame axes. 

th 

FP(I,J) — is the I component of the external force acting upon 
particle J in the appendage (I = 1,2,3; J = 1,2,.. ,,N). 

The external forces for each of the particles comprising the appendage 
are resolved along body axes. 

Subroutines 

Subroutine INITL reads in all program input data and performs 
consistency checks. Selected input data is echo printed. The eigenvalues 


* 

denotes the nirnber zero. 
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and eigenvectors of the standard syxcunetric eigenvalue problem grven by 
Eq. (5-2) are computed via a call to IMSL subroutine EIGRS. The eigen- 
vectors are then transformed to those corresponding to Eg. (5-1). All 
time-invariant terms of the generalized mass matrix of Figure 3 are 
computed. Finally, the initial conditions on the particle displacements, 
modal coordinates, and the respective time derivatives are set. 

Subroutine GNMASS computes and assembles the deformation depend 
ent generalized mass matrix of Figure 3. 

Subroutine EXTF computes and assembles the generalized force 
vector of Eq. (8-2) . 

Subroutine NLKT computes and assembles the vector of nonlinear 
kinematic terms of Eq. (8-3) . 


subroutine SOLVE computes the transformation matrices given by 
Eq. (7-1) and (7-2). The set of simultaneous equations given by Eq. (8-1) 
are solved via a call to IMSL subroutine LEQTlP. The state vector 
Eq (8-4) is assembled and its time derivative, (8-5), evaluated. 

The value of the state vector is advanced one time step via a call to 

subroutine ODESLV. 

subroutine ODESLV integrates the state equation, Eq. (8-5), using 
the Adams method with third order differences. 

subroutine PEim is executed only at print-time intervale speci- 
fied in the input (see below) . When called, the subroutine prints the 
time, force, and torque on the rigid body, applied forces on the par- 
ticles and all the variables of the vector given in Eq. (9 1). 


Entry point GHEE in subroutine PRIET stores selected variables 
for plotting at a specified time Interval (see naselist items DTG and 

IPLOT below) . 
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Program Input Data 


Pj^ogram input data is r©ad in during exscution of subroutino INITL. 
Input is achieved through four READ- NAMELIST combinations and a single 
unformatted READ of the stiffness matrix. It is worth noting that while 
the code given in Appendix A requires the stiffness matrix (described in 
Section 3.1) and from this and the appendage mass matrix (assembled in- 
ternally) computes the constrained appendage eigenvalues and eigenvectors, 
it could be modified to read in the appropriate eigenvalues/eigenvectors 
directly. The four NAMELIST inputs are defined below, and their use 
illustrated in Appendix B. 

(1) NAMELIST/INPUT /M(?, N, MASS, RM, 10, S, NT; contains all 
mass and geometry data as well as the number of modes to 
be retained. 


M0 

N 


MASS 

RM 


10 


S 


NT 


mass of rigid body (real) 
number of particles (integer) 

masses of particles 1 through N (real N x 1 array) 

position vectors of particles 1 through N prior 
to deformation, expressed in body frame (real 
3 X N array) 

inertia matrix of the rigid body with respect to 
a frame located at the rigid body mass center 
with axes parallel to body frame (real 3x3 
array) 

position vector from body frame origin to mass 
center of rigid body expressed in body frame (real 
3 x 1 array) 

number of modes to be retained; modes 1 through NT 
are used (integer) 
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(2) NAMELIST/KIN/UVW, OMEGA, R, THETA: contains initial con- 

ditions for kinematic variables. 

UVW = initial velocity vector of body frame orxgrn, 

expressed in body frame coordinates (real 3 x 1 

array) 


OMEGA = initial angular velocity vector of body frame 
with respect to inertial frame, components ex 
pressed in body frame (real 3 x l array) 



R 

= initial inertial position vector of body frame 
origin, components expressed in inertial frame 



(real 3 x 1 array) 


THETA 

_ initial 1-2-3 Euler angles of body frame with 
respect to inertial frame (real 3 x 1 array) 

(3) 

NAMELIST /RUN /DT, TSTOP, DTP, DTG: contains numerical 

integration parameters and print and plot time intervals. 


DT 

= integration time step in seconds (real) 


TSTOP 

= integration termination time in seconds (real) 


DTP 

= print output time interval in seconds; output 
printed every DTP seconds (real) 


DTG 

= plot output time interval in seconds: selected 
variables plotted every DTG seconds (real) 

(4) 

NAMELIST/PLOT/IPLOT: specifies which elements of vector 

in Eq. (9-1) are to be plotted via line printer. 


IPLOT = integer array with the integers corresponding 
to those elements of the vector in Eg. (9-1) 
that are to be plotted versus time (see sample 


use in Appendix B) . 


i 

f 

I 
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APPENDIX A 


FORTRAN PROGRAM LISTING 


LEVEL Z.3.0 (JUKE 78) OS/3SO FORTRAN H EXTCMOED DATE 82.271/12.^3.05 

REQUESTED OPTIONS-* N0C3J .TERN » ,K3XREF , »NONAP. » ,NAME( MAIN) , AO (NONE ) ,OPT( 0 , , FLAG( I ) .SIZE ( 30^iK ) , LC( 60 ) , 

OPTIONS IN effect: NAME(MAIN) NOQPTIMIZE IIMECOUNTC 60 ) SIZE(0360K) A'JTCDOUMOKE ) 

SOURCE EBCDIC NOLIST NOOECK NOOBJECT NONAP NOFORIIAT C03TMT NOXREF NOALC NOANSF TERM IDM 


ISN 0002 
ISN 0003 
ISN 0006 
ISN 0005 


ISN 0006 
ISN 0007 
ISN 0008 


ISN 0009 


ISN 0010 


ISN 0011 


ISN 0012 
ISN 0013 
ISN 0016 


C 

C 

C 

10 

c 

c 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


OOOODIOO 

•THIS PROGRAM SOLVES THE EQUATIONS OF MOTION OF A VEHICLE 00003300 

• COHSISTir.'S OF A RIGID BASE WITH AH ATTACHED FLEXIBLE APPENDAGE. 00000600 

•THE APPEK’BAGE IS IDEALIZED AS A COLLECTION OF PARTICLES CO tNECTED OOOOOEOO 

•or MASSLESS ELASTIC STRUCTURE. IN ADDITION TO EXTERNAL FORCES ACTINGOOOOO LOO 
•UPON EACH OF THE PARTICLES » A FORCE AIJD TORQUE ARE ACCOKfrODATED AT 00000700 
•THE GRAPPLE FIXTURE CORRESPCirDIMG TO THE ORIGIN OF BODY FRAME.. 00000800 

•njRITTEN BY JOEL 5TORCH A STEPHEr( GATES C.S.O.L. BASED UPON 00000900 

• C.S.O.L. REPORT CR15G2 SEPT. 1982) 00001000 

00001200 

00001300 

note: arrays ARE DIMENSIONED TO ACCOMODATE A MAXIMUM OF 50 PARTICLES00001^«00 

00001500 

IMPLICIT REAL"8( A-H,0-Z) 

OIMENjICN F0( 3),TAU0( 3)»FP( 3,50) 

COMKON /TIME/ DT.TSTOP ►DTP ,OTG 
T = 0.0 

INPUT PROGRAM DATA AND CALCULATE ALL TIME INVARIANT PARAMETERS. 


CALL INITL 

TP=OTP 

TG=OTG 

CALCULATE GENERALIZED MASS MATRIX 
CALL GhaiASS 

INPUT VALUES REQUIRED FOR *EXTF' 

(ALL VECTORS EXPRESSED IN BODY FRAME) 

FO - EXTERNAL FORCE ON NO (AT ORIGIN OF BODY FRAME) 

TAUO - EXTERNAL TORQUE AT LOCATION OF BODY FRAME ORIGIN. 
FP - VECTOR OF EXTERNAL FORCES ON PARTICLES I.2»...»N. 

CALL EXTF(F0,TAU0»FP) 

CALCULATE NON-LINEAR KINEMATIC TERNS 


CALL NLKT 

INTEGRATE EQUATIONS OF MOTION 

CALL SOLVE 
T=T*DT 

IF(T .LT. TP) GO TO 15 


OOOOlCiOO 

OC001700 

03001600 

00001900 

0C002000 

00002100 

OC002200 

00002300 

00002600 

00002500 

00002600 

00002700 

00002880 

00002900 

00003000 

00003100 

00003200 

00003300 

00003600 

00003500 

00003600 

00006300 

OOU06600 

OQ006SOO 

00006600 

00006700 

00006800 

00006900 

00005000 

00005100 

00005200 

00005300 


PAGE 1 


FLAG( I ) 
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ruiN 


OS/360 FORTRAN H EXTENDED 


DATE 82.2^6/12.21.^^ 


PAGE 2 


LEVEL 2.3.0 CJUNE 


78) 


ISN 0026 
I5N 0027 
ISN 0028 
ISN 0030 
ISN 0031 
ISN 0032 
ISN 003^ 
ISN 0035 
ISN 0036 


CALL PRINTlT,F0»TAU0,FP) 
TP=TP+0TP 

15 IF(T .LT. TG) SO TO 20 
CALL GRAFCT) 

T6=TG+DTG 

20 IFIT .LT. TSTOP) GO TO 10 
CALL PLOT 
STOP 
EW) 


00005400 

00005500 

00005600 

00005700 

00005800 

00005900 

00006000 

00006100 

00006200 


:::: : ™ - 

PPOGPAH SIZE = 2002 » SUBPROGRAM NAME = MAIN 

.STATISTICS. SOURCE STATEMENTS * 35. PROGRAM SIZE 


I) 


•STATISTICS* NO DIAGNOSTICS GENERATED 
END OF COMPILATION ****** 


280K BYTES OF CORE NOT USED 
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DATE dZ.2<t6/12.21.4S 


PAGE 1 


LEVEL 2.3.0 UUNE 78) 03/360 FORTRAN H EXTENDED 

REQUESTED OPTIONS: NOOBJ.TERM. »NOXREF » »NOttAP» p pNAME(MAIN) p AOCNONE ) pOPTt 0 ) p p pFLAGC I J pSIZEt 38^K ) p LCt 60 ) p 

OPTIONS IN effect: NAME (MAIN) NOOPTIMIZE LINECOUNTl 60 ) SIZE(0384K) AUTOOBLtNONE ) 

SOURCE EBCDIC NOLIST NOOECK NOOBJECT NOMAP HOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM FLAG(I) 


ISN 0002 SUBROUTINE SKEW(VpA) 00006300 

Q 00006400 

C THIS SUBROUTINE CREATES THE SKEW SYMMETRIC MATRIX CORRESPONDING 00006500 

C TO THE VECTOR "V". 00006600 

Q 00006700 

ISN 0003 REAL»6 VpA ' 00006600 

ISN 0004 DIMENSION V(3)pA(3p3) 00006900 

ISN 0005 A(1p1)=0.0 . 00007000 

ISN 0006 A(1p2)=-VC3) 00007100 

ISN 0007 A(1,3)=V(2) 00007200 

ISN 0008 A(2,X)=V(3) 00007300 

ISN 0009 A(2p2)=0.0 00007400 

ISN 0010 A(2p3)=-V(1) 00007500 

ISN 0011 A(3,1)=-V(2) 00007600 

ISN 0012 A(3,2)=V(1) 00007700 

ISN 0013 A(3p3)=0.0 00007800 

ISN 0014 RETURN 00007900 

ISN 0015 END 00008000 

•OPTIONS IN EFFECT»NAME(MAIN) NOOPTIMIZE LINECOUNT(60 ) SIZE(0384K) AUTOOBL(NONE ) 

•OPTIONS IN EFFECT*SOURCE EBCDIC HOLIST NOOECK NOOBJECT NOMAP NOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM FUGCI) 
•STATISTICS^ SOURCE STATEMENTS = 14 p PROGRAM SIZE = 388, SUBPROGRAM NAME = SKEW 

•STATISTICS* NO DIAGNOSTICS GENERATED 


i m i n i w w END OF COMPILATION ••<hhi* 280K BYTES OF CORE NOT USED 
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LEVEL 2.3.0 (JUNE 
requested options: 


78) 

NOOBJ.TERM» »NOXREFi tNOfUPi 


OS/360 FORTRAN H EXTENDED 

, »NANE( MAIN ) ,A0( NONE ) ,OPT( 0 


date 82. 246/12. 21. A6 
, , FUGl I ) . SIZE! 384K ) , LCC 60 ) » 


PAGE 


OPTIONS IN EFFECT 


; NAHEmAIN) NOOPTINIZE 
SOURCE EBCDIC NOLIST 


l 1NFC0UNT( 60 ) SIZE(0384K) AUTOOBLC NONE ) 

.c^T KinMAP noformat gostwt noxref 






TERN IBN FLAG(I) 


ISN 0002 


ISN 0003 
ISN 0004 
ISN OOOS 
ISN 0006 
ISN 0007 
ISN 0008 


SUBROUTINE CROSSCA»B.C) 

C fflis SUBROUTINE CALCULATES THE VECTOR CROSS PRODUCT 
C A X B =C 

c 

REALK8 A(3J»B(3)»C(3) 

C(1)=A(2)*B(3)-A(3)*»BC2) 

C(2)=A(3)*‘Bll)-A(l)**B(3) 

C( 3 )=A( 1 1»B( 2 )-At 2 )«B( 1 ) 

RETURN 

END 


00008100 

00008200 

00008300 

00008400 

00006500 

00008600 

00008700 

00008800 

00008900 

00009000 

00009100 


.OPTIONS IN EFPECT.KAHE.HAIK. NOOPTIHIZE LINECOUHT. 60 . SIZE.0S.6K, AUTOOBL.HONE. 

Ti:: IN EFFECT.SOURCE EBCDIC .LIST .DECK ^.ECT .XR.P .FORHAT OOSTHT .«XREF .ALC NOANSF TERH IBH FLAO.I. 
.STATISTICS. SOURCE STATEHENTS = 7. PROORAH SIZE = SUBPROORAH NAHE = CROSS 


»5jj^TXSTICS» NO DIAGNOSTICS GENERATED 
•****» END OF CONPIUTION ****** 


280K BYTES OF CORE NOT USED 


1 
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LEVEL 2.3,0 (JUNE 78) 


03/360 FORTRAN H EXTENDED 


DATE 82.2^6/12.21.67 


PAGE 


REQUESTED OPTIONS: N008J,TERM, ,NOXREF , »NOtlAP, , ,NAME( MAIN ) , A0( NONE ) »OPT( 0 ) » » »FLAG( I ) »SIZE( 386K ) , LC(60 ) , 

OPTIONS IN effect: NAME(MAIN) NOOPTIMIZE LINEC0UNT( 60 ) SIZE(0386K) AUTOOBLC NONE ) 

SOURCE EBCDIC NOLIST NODECK NOOBJECT NOMAP NOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM FLAG(I) 


ISN 0002 


ISN 0003 
ISN 0006 
ISN 0005 
ISN 0006 


ISN 0007 
ISN 0008 
ISN 0009 
ISN 0010 
ISN 0011 
ISN 0012 
ISN 0013 
ISN 0016 
ISN 0015 
ISN 0016 
ISN 0017 
ISN 0016 
ISN 0019 
ISN 0020 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE INITL 

THIS SUBROUTIt(E READS IN PROGRAM DATA AND CALCUUTES ALL TIME 
INVARIANT PARAMETERS. 

IMPLICIT REAL«8(A-H»0-Z) 

REAL«e MASS tMO , K ,MC0 , INERTl ,MRM , 10 >M0S 
REAL«6 PLTDAT 

DIMENSION MASS( 50 ) » RM( 3 , 50 ) > K( 150 , 150 ) >MC0 ( 3 ) > INERTl ( 3 * 3 )> 

1 MRM(3»50)»CMAT(3,3).A12(3»3),A23(3»150),A(156>156)>Q(3>50). 

1 O00T(3t50)tUVU(3),OMEGA(3),R(3)»THETA(3).ia(3>3).S(3)>M0S(3)> 

2 AV( 11325 ) > FREQ( 150 ) » EV( 150 * 150 ) »HK( 150 ) > IPLOT( 62 ) * PLTDATC 100,20), 

3 t4KM(3,150),ETA(150),ETAD(150) 

COMMON /CONST/ TM ,MASS,MC0 ,MRM,CMAT,N,N3,N3P6 ,KT,NTP6 ,NO 

COMMON /AMAT/ A, A12, INERTl, A23 

COMMON /STATE/ R , THETA, Q.UVW, OMEGA, QDOT 

COMMON /GEOM/ RM 

COMMON /TIME/ OT,TSTOP,OTP,DTG 

COMMON /RIGID/ M0S,I0,S 

COMMON /PLOTT/ P LTD AT, IP LOT, NP 

COMMON /MOOCO/ ETA,ETAD 

COMMON /MOOES/ EV,FREQ 

EQUIVALENCE ( K( 1 , 1 ) , EV( 1 , 1 ) ) 

NAMELIST /INPUT/ MO,N,MASS,RM,IO,S,MT 
NAMELIST /KIN/ UVW, OMEGA, R, THETA 
NAMELIST /RUN/ OT,TSTOP,DTP,OTG 
NAMELIST /PLOT/ IPLOT 

DESCRIPTION OF NAMELIST VARIABLES 

••MO" IS THE MASS OF THE RIGID BASE 

••N’* IS THE NUMBER OF PARTICLES THAT COMPRISE THE FLEXIBLE APPENDAGE 
"MASS** CONTAINS THE MASSES OF PARTICLES 1 TO N 

••RM" CONTAINS THE POSITION VECTORS OF PARTICLES 1 THRU N IN THE 
UNOEFORMED STATE EXPRESSED IN THE BOOT FRAME. 

"10" IS THE INERTIA MATRIX OF THE RIGID BASE WITH RESPECT TO A 
FRAME LOCATED AT THE MASS CENTER OF THE BASE AND PARALLEL 
TO THE BODY FIXED AXIS SYSTEM. 

"S" IS THE VECTOR FROM THE BODY FRAME ORIGIN TO THE MASS CENTER 
OF THE RIGID BODY. 

••NT" NUMBER OF RETAINED MOOES IN APPENDAGE VIBRATION 

••UVW" IS THE INITIAL VELOCITY OF THE BODY FRAME ORIGIN EXPRESSED 
IN BODY COORDINATES. 

"OMEGA'* IS THE INITIAL ANGULAR VELOCITY OF THE BODY FRAME EXPRESSED 
IN BODY COORDINATES. 


00009200 

00009300 

00009600 

00009500 

00009600 

00009700 

00009800 

00009900 

00010000 

00010100 

00010200 

00010300 

00010600 

00010500 

00010600 

00010700 

00010800 

00010900 

00011000 

00011100 

00011200 

00011300 

00011600 

00011500 

00011600 

00011700 

00011800 

00011900 

00012000 

00012100 

00012200 

00012300 

00012600 

00012500 

00012600 

00012700 

00012800 

00012900 

00013000 

00013100 

00013200 

00013300 

00013600 

00013500 

00013600 

00013700 

00013800 

00013900 

00016000 

00016100 

00016200 

00016300 

00016600 
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LEVEL 2.3.0 (JUNE 78) 


INITL 


OS/360 FORTT?AN H EXTENDED 


ISN 0039 
ISN 00^0 
ISN 0041 
ISN 0042 
ISN 0043 
ISN 0045 
ISN 0046 
ISN 0048 
ISN 0049 
ISN 0050 


ISN 0051 
ISN 0052 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


15 THE INITIAL INERTIAL POSITION VECTOR OF THE BOOT FRAME 
ORIGIN 


ISN 0021 
ISN 0022 
ISN 0024 
ISN 0025 
ISN 0026 

8 

ISN 0027 
ISN 0028 

10 

ISN 0029 

100 

ISN 0030 

101 

ISN 0031 
ISN 0032 
ISN 0033 
ISN 0034 
ISN 0035 

103 

ISN 0036 
ISN 0037 

104 

ISN 0038 



■•THETA” IS THE INITIAL SET OF ATTITUDE ANGLES FOR THE BOOT FRAME 

'*OT** IS THE INTEGRATION TIME STEP 

••TSTOP" IS THE TERMINAL TIME FOR THE SIMULATION 

"OTP” IS THE TIME INTERVAL BETWEEN PRINTOUTS 

"DTG” IS THE TIME INTERVAL BETWEEN PLOTTED POINTS 

••IPLOr* IS AN ARRAY INDICATING VARIABLES TO 

numbering CORRESPONDS TO LOCATION IN STATE VECTOR. 
POINTS ARE PLOTTEO EVERY “OTG- SECONDS. 

READ(5, INPUT) 

IFINT .LE. 3*N) GO TO 8 
WRITE(6iI12) NT,N 
STOP 

WRITEC6»100) MO 

DO 10 I=lfN . . , ,5 

WRITE(6,101) IpMASSm*(RM(JRl).J»lt3) ✓/ IH 

FORMAT(1H1.20X.*MASS OF RIGID BASE=SE12.4* SLl^ »//»lH 

I 'PARTICLESTEl.'MASStSLUGSJST^l.’POSITIWCFT.) »/) 

FORMATCIH ,T12»I3.T21»F7.2fT34,3(F7.2.2Xn 
WRITE(6il06 ) I ( I0( It J)» J-l»3)»I“l»3i 
WRITE(6,107) SiNT 
REA0(5.KIN) 

?SrMAT?1ho!5X^ INITIAL 

1 -INITIAL ANGULAR VELOCITT=‘ .3F7.2. ‘ OEG/SEC > 

1 


DEG* ) 


'Z ISN 0053 


108 

109 


110 

105 


112 

106 


107 


INITIAL ATTITLrt)E = * »3F8.3» ' 

READ<5»RUN) 

WRITE(6>105) OT»TSTOP»0TP»DTG 
REA0(5»PL0T) 

11*0 

DO 108 1=1*42 

IFdPLOTd) .EQ. 0) 60 TO 109 
II*II^1 

IFdl .EQ. 0) 60 TO 111 
WRITEt6*llO) tlPLOTd). 1=1*1^ 

1 ElZ.'t.' SEC .2X. -PRINT INTERVAL^- .E12-4. ‘ SEC . PLOT INTERVAL . 

Sci^TUHO^IxJlZ,- MOOES RE0UESTE0-.2X.I3,- PARTICLES IN IME*-- ) 

FORlJiTl lHSlnS^-iNERTIA MATRIX OF RIGID BOOYCSLUG FT«2)- . 

1 //»<T1Si 3E13.5) ) ^ 3XtI2*' CONSTRINED APPENDAGE MOOE00019900 

F0RMAT(1H0»15X» *S- ,3E13.5* FT 00020000 

IS RETAINED* ) 00020100 

CHANGE ANGULAR VELOCITY A ATTrOOE TO RADIAN MEASURE 00020200 


DATE 82.246/12.21.47 


00014500 
00014600 
00014700 
00014800 
00014900 
00015000 
00015100 
00015200 
00015300 
00015400 
00015500 
00015600 
00015700 
00015800 
00015900 
00016000 
00016100 
00016200 
00016300 
00016400 
00016500 
00016600 
00016700 
00016800 
00016900 
00017000 
00017100 
00017200 
00017300 
00017400 
00017500 
00017600 
00017700 
00017800 
00017900 
00018000 
00018100 
00018200 
00018300 
00018400 
00018500 
00018600 
00018700 
00018800 
00018900 
00019000 
00019100 
00019200 
00019300 
00019400 
00019500 
00019600 
00019700 
00019800 


.no. 


PAGE 2 


44 


LEVEL 2.3,0 (JUNE 781 


INITL 


OS/360 FORTRAN H EXTENDED 


DATE 82 , 2 ^ 6 / 12 . 21 . 


PAGE 3 


C 


ISN 0054 

111 

0TR=DATAN( 1 ,0D0 )/45. 

ISN 0055 


DO 15 I=l>3 

ISN 0056 


OMEGA( I )=DTR*OMEGA( I ) 

ISN 0057 

15 

THETA ( I )=0TR*THETA( I ) 


C 

c 

TM - TOTAL BODY MASS 

ISN 0058 

u 

TM=M0 

ISN 0059 


DO 20 1=1 iN 

ISN 0060 

20 

TM=TM^MASS(I) 

ISN 0061 


DO 30 J=l»3 

ISN 0062 


M0S( J)=M0»S( J) 

ISN 0063 

30 

MC0( Jl=0.0 

ISN 0064 


DO 40 I=liN 

ISN 0065 


DO 45 J=lf3 

ISN 0066 


MRM( J , I )=MASS( I )«RM( J , 1 ) 

ISN 0067 

45 

MC0( J )=MC0( J )+MRM( J »I ) 

ISN 0068 

40 

CONTINUE 

ISN 0069 


DO 42 I=l>3 

ISN 0070 

42 

MC0(I)=MC0(I)*M0S(I) ^ 

ISN 0071 


N3=3»N 

ISN 0072 


N3P6=N3«6 

ISN 0073 


NTP6=NT^6 

ISN 0074 


NO=2«NTP6 

ISN 0075 


NP=0 


c 

C 

READ IN STIFFNESS MATRIX 

ISN 0076 

c 

REA0(8) NDOF,(CK(I,JliJ=l»NDOF),I=l»NDOF) 

ISM 0077 


IF(NDOF .EQ. N3) GO TO 400 

ISN 0079 


WRITE(6.102) N»NDOF 

ISN 0080 


STOP 

ISN 0081 

400 

CONTINUE 

ISN 0062 

102 

FORMATdHO.lOX.* INCONSISTENT DATA* .2X.I3» * PARTICLES* .2X. 13. 



1 • DEGREES OF FREEDOM IN STIFFNESS MATRIX* ) 


C 

c 

GET CONSTRAINED FREQUENCIES AND MODE SHAPES OF APPENDAGE 

ISN 0083 

c 

. L=1 

ISN 0084 


DO 300 J=1»N3 

ISN 0085 


LC=l+J/3 

ISN 0086 


IF( (J-3»»(J/3n .EQ. 0) LC=LC-1 

ISN 0068 


DO 300 1=1. J 

ISN 0089 


LR=l+I/3 

ISN 0090 


IFC (I-3*( 1/3)1 ,EQ. 0) LR=LR-1 

ISN 0092 


AV( L )=K( I p J )/DSQRT( MASS! LR )*MASS( LC ) ) 

ISN 0093 


L=L*1 

ISN 0094 

300 

CONTINUE 

ISN 0095 


CALL E IGRS ( AV . N3 , 1 , FREQ , EV , 150 , . lER ) 

ISN 0096 


IFdER .EQ, 0) GO TO 310 

ISN 0098 


HRITE(6,301) lER 

ISN 0099 

301 

FORMATdHO.lOX, 'ERROR FROM IMSL ROUTINE "EIGRS** ERROR C00E=',I4) 

ISN 0100 


STOP 


C 

C TRANSFORM EIGENVECTORS 
C 

ISN 0101 310 DO 311 L=liN 


00020300 

00020<t00 

00020500 

00020600 

00020700 

00020800 

00020900 

00021000 

00021100 

00021200 

00021300 

00021400 

00021500 

00021600 

00021700 

00021800 

00021900 

00022000 

00022100 

00022200 

00022300 

00022400 

00022500 

00022600 

00022700 

00022800 

00022900 

00023000 

00023100 

00023200 

00023300 

00023400 

00023500 

00023600 

00023700 

00023800 

00023900 

00024000 

00024100 

00024200 

00024300 

00024400 

00024500 

00024600 

00024700 

00024800 

00024900 

00025000 

00025100 

00025200 

00025300 

00025400 

00025500 

00025600 

00025700 

00025800 

00025900 

00026000 
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LEVEL 2.3.0 (JUNE 781 


OS/360 FORTRAN H EXTEHOEO 


ISN 0102 
ISH 0103 
ISN 01Q<* 
ISN 0105 
ISN 0106 
ISN 0107 
ISN 0108 
ISH 0109 
ISN 0110 

ISN 0111 


ISN 0112 
ISN 0113 
ISN 011^ 
ISN 0115 
ISN 0116 
ISN 0117 
ISN 0118 
ISN 0119 
ISN 0120 
ISN 0121 
ISN 0122 
ISN 0123 
ISN 0129 
ISN 0125 
ISN 0126 
ISN 0127 
ISN 0128 
ISN 0129 
ISN 0130 
ISN 0131 
ISN 0132 
ISN 0133 
ISN 0139 
ISN 0135 
ISN 0136 
ISN 0137 
ISN 0136 
ISN 0139 
ISN 0190 
ISN 0191 
- ISN 0192 
ISN 0199 
ISN 0195 
ISN 0196 


C1=0SQRT(MASS(D) 

00 311 1=1.3 
IR=3»(L-1)*I 
DO 311 J=1.N3 
EV(IR.J)=EV(IR,J)/C1 
DO 320 1=1. N3 

FHZ=0SQRT( FREQC I ) )/6 . 283185 

1 'MODE SHAPE! • , ( 9(S12 .9.1X) ) ) 

CONTINUE 

Tue qrmTTNF "EIGRS*' RETURNS AN ORTNONORMAL SET OF EIGENVECTC3RS. 

SlScrSrASSUriE in the DERmTICN TH*T TOE 
ElfiENVECTORSI OF TOE ORIGINAL GENERALIZED EIGENVALUE PR06LEM) 
«1 totoS? NITO RESPECT TO HASS« STIFFNESS ) AND HORNALIZEO 
HITH RESPECT TO MASS. 

CALCUUTE TIME INVARIANT PART OF INERTIA MATRIX "INERTl” 
and "CMAr* 

00 50 1=1.3 

DO 50 J=I»3 

INERT1(I.J)=0.0 

CMAT(I.J)=0.0 

CONTINUE 

DO 60 1=1 .N 

XS=RMC1,I)**2 

YS=RM(2.I)»**2 

ZS=RM(3.I)***Z 

INERTm,ll=IHERTm.ll.MASS(It»«TS»ZS| 

TNERTlt 1 ,2 )=INERT1(1 pZ )-HASS( II*RM(lf I )*RM( Ztll 

SEm(i:r.=TOiRTi(i:3i-«Ass(i.w^^ 

INERTK 2f2)=INERTl(2.2)^MASS(l)»(XS>ZS) 

INERTl (2.3 ) =INERT1 (2.3) -MASS( I )»RM( 2 . I )»RM( 3.1) 

INERTK 3 . 3 )=INERT1( 3,3) 4MASS( I )*( XS* YS ) 

CMAT( 1 , 1 )=CMAT( 1.1) ♦MASSC I )*»XS 

CHAT* 2 . 2 ) =CMAT( 2 , 2 ) ♦MASS! I )»YS 

CMATC 3 . 3 )=CMAT( 3 , 3 )4MASS( I )«ZS 
) CONTINUE 

CMAT( 1 . 2 ) s-INERTK 1.2) 

CMAT(1,3)=-INERTI(1.3) 

INERTl (1,2 )=INERT1 (1.2)410(1,2) -M0»S( 1 2 ) 

INERTK 1,3 ) = INERT1( 1,3 )4l0( 1.3 )-M0*»S( I )»S( 3) 

INERTl(2,2)=INERTl(2,2)4l0(2.2)4M0tKS(l)^24S(3)**2) 

INERTK 2.3) = INERT1( 2.3)4l0( 2.3)-M0»S( 2)*S( 3) 

wIr^kL3) = IHERTK3,3)4I0(3,3)4M0MS(1)»^^^ 

DO 70 1=1.3 
DO 70 J=1.3 
IFd .LE. J I 60 TO 70 
INERTK I . J ) = INERTK J. I ) 

CMAT(I.J)=CMAT(J,I) 

0 CONTINUE 


)4M0»( S( 2 )«»24S( 3 )**2 ) 
)-M0»S( 1)*S( 2) 
)-M0*»S(I)»S(3) 

)4M0»( S( 1 )**24S( 3 )*»2 ) 

,)-nO*S(2)*S(3) 

1 )4M0*»( S( 1 )*»»24S( 2 )«H*2 ) 


CREATE TIME INVARIANT PORTIONS OF 6EHEHALIZE0 MASS MATRIX "A" 


DATE 82.296/12.21.97 

00026100 

00026200 

00026300 

00026900 

00026500 

00026600 

00026700 

00026800 

00026900 

00027000 

00027100 

00027200 

00027300 

00027900 

00027500 

00027600 

00027700 

00027800 

00027900 

00028000 

00028100 

00028200 

00028300 

00028900 

00028500 

00028600 

00028700 

00028800 

00026900 

00029000 

00029100 

00029200 

00029300 

00029900 

00029500 

00029600 

00029700 

00029800 

00029900 

00030000 

00030100 

00030200 

00030300 

00030900 

00030500 

00030600 

00030700 

00030600 

00030900 

00031000 

00031100 

00031200 

00031300 

00031900 

00031500 

00031600 

00031700 

00031600 
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PAGE 



c 


00031900 


c 

(1,2) PARTITION "A12" 

00032000 


c 


00032100 

ISN 0167 


CALL SKEN(riC0,A12) 

00032200 

ISN 0168 


DO 80 1=1,3 

00032300 

ISN 0169 


DO 60 J=l,3 

00032600 

ISN 0150 


IF(I .EQ. J) GO TO 80 

00032500 

ISN 0152 


A12(I,J)=-A12(I,J) 

00032600 

ISN 0153 

80 

CONTINUE 

00032700 


c 


00032800 


c 

(2,3) PARTITION "A23" 

00032900 


c 


00033000 

ISN 0156 


DO 90 1=1, N 

00033100 

ISN 0155 


L=3«I-2 

00033200 

ISN 0156 

90 

CALL SKEW(NRn(l.I),MKn(l,L)) 

00033300 

ISN 0157 


DO 92 1=1,3 

00033600 

ISN 0158 


DO 92 J=1,NT 

00033500 

ISN 0159 


A23(I,J)=0.0 

00033600 

ISN 0160 


00 96 L=1,N3 

00033700 

ISN 0161 

96 

A23(I,J)=A23(I,J)+WKmi,L)**EVCL,J) 

00033800 

ISN 0162 

92 

CONTINUE 

00033900 


c 


00036000 


c 

STORE CONSTANT PARTITIONS OF “A" 

00036100 


c 


00036200 

ISN 0163 


00 200 I=1,HTP6 

00036300 

ISN 0166 


00 200 J=1,NTP6 

00036600 

ISN 0165 


A(I,J)=0.0 

00036500 

ISN 0166 

200 

CONTINUE 

00036600 


c 


00036700 


c 

CREATE (1,1) PARTITION 

00036800 


0 


00036900 

ISN 0167 


DO 210 1=1,3 

00035000 

ISN 0168 


DO 210 J=I,3 

00035100 

ISN 0169 


IF(I .EQ. J) A(I,J)='ni 

00035200 

ISN 0171 

210 

CONTINUE 

00035300 


c 


00035600 


c 

CREATE (1,3) PARTITION 

00035500 


0 


00035600 

ISN 0172 


DO 215 1=1,3 , 

00035700 

ISN 0173 


DO 215 J=1,N3 

00035800 

ISN 0176 

215 

NXM(I,J)=0.0 

00035900 

ISN 0175 


00 220 L=I,N 

00036000 

ISN 0176 


JS=3*L-2 

00036100 

ISN 0177 


DO 221 1=1,3 

00036200 

ISN 0178 


MKmi,JS)=HASS(L) 

00036300 

ISN 0179 


JS=JS+1 

00036600 

ISN 0180 

221 

CONTINUE 

00036500 

ISN 0161 

220 

CONTINUE 

00036600 

ISN 0182 


DO 283 1=1,3 

00036700 

ISN 0183 


DO 283 J=1,NT 

00036800 

ISN 0186 


JP6=J^6 

00036900 

ISN 0185 


A(I,JP6)=0.0 

00037000 

ISN 0186 


DO 281 L=1,N3 

00037100 

ISN 0187 

281 

A(I,JP6)=A(I,JP6 )♦«<«( I, L)*EV( L,J) 

00037200 

ISN 0188 

263 

COt(TINUE 

00037300 


0 


00037600 


C 

CREATE (3,3) PARTITION 

00037500 


c 


00037600 
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LEVEL 2»3.0 lJUNE 76) 
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ISN 0169 
ISN 0190 


ISN 0191 
ISN 0192 
ISN 0193 
ISN 0196 
ISN 0195 
ISN 0196 
ISN 0197 
ISN 0196 
ISN 0199 
ISN 0200 
ISN 0201 


DO 230 L=1»NT 
230 A( L*6*L+6 )=l. 

c SET INITIAL DEFOPrUTIOH AND RATE TO ZERO 

c 

DO 250 1^1 »N 
DO 250 J=l»3 
qu>X)=o.o 

QOOTU,I)=0.0 
250 CONTimiE 

00 293 1=1 iNT 
ETA(I)=0.0 
ETAD(I)=0.0 
293 CONTIWE 
RETURN 
END 


DATE 62-266/12.21.67 

00037700 

00037600 

00037900 

00036000 

00036100 

00036200 

00036300 

00036600 

00036500 

00036600 

00038700 

00036800 

00036900 

00039000 

0003910Q 

00039200 


ISN OZOi 

PROGRAM SIZE * 1030Z6, SUBPROGRAM NAME = INITL 

.STATISTICS. SOURCE STATEMENTS = 200, PROGRAM SIZE 

»«t.t1STICS. no DIAGNOSTICS GENERATED 
.STATISTICS 

,«.« END OF COMPILATION .««* 
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DATE 02.2<»6/12.21.50 


PAGE I 


REQUESTED OPTIONS: NOOBJ »TERM, »NOXREF 


>NOMAP» I ,NAME(MAIN) lAD(NONE) lOPTt 0 ) » » tFLA6( I ) >SIZ£( 36«iK)»LC(60 ) » 


OPTIONS IN effect: 


NAME(rUIN) NOOPTIMIZE LINECOUNT( 60 ) SIZE(03SAK) AUTODBL( NONE ) 
SOURCE EBCDIC NO LIST NODECK NOOBJECT NONAP NOFORHAT GOSTm NOXREF 


NOALC NOANSF TERM IBM FLAG(I) 


ISN 0002 


ISN 0003 
ISN 0004 
ISN 0005 


ISN 0006 
ISN 0007 
ISN 0008 
ISN 0009 
ISN 0010 
ISN 0011 
ISN 0012 
ISN 0013 
ISN 0014 
ISN 0015 
ISN 0016 
ISN 0017 


ISN 0018 
ISN 0019 
ISN 0020 
ISN 0021 
ISN 0022 
ISN 0023 
ISN 0024 
ISN 0025 
ISN 0026 
ISN 0027 
ISN 0028 
ISN 0029 
ISN 0030 
ISN 0031 


ISN 0032 
ISN 0033 
ISN 0034 
ISN 0035 
ISN 0036 


ISN 0037 
ISN 0038 
ISN 0039 


SUBROUTINE GNMASS 

C THIS SUBROUTINE CALCULATES THE GENERALIZED MASS MATRIX '*A" 

C 

IMPLICIT REAL»8(A-H,0-Z) 

REALMS MA5S,MC0,MRM.INERTI»MQ»MC1,INERT2,INERT3 
DIMENSION MASS(50 ).MC0( 3) »MRM( 3»50)fCMAT(3»3) rA( 156*156 )> 

1 A12(3,3),INERT1(3,3)»A23(3»1501,Q(3,50)»MQ(3*50),MC1(3). 

2 KN(3,3),INERT2(3*3),A23q(3,150),INERT3(3,3),QOOT(3.50),UVWC3)* 

3 0MEGA(3),R(3>*THETA(3) 

COMMON /CONST/ TM»MASS*MC0 *MRM,CMAT*N*N3*N3P6 *NT,NTP6 *NO 

COMMON /AMAT/ A* A12 » INERTl » A23 

COMMON /STATE/ R * THETA ,q,UVW* OMEGA, QOOT 

COMMON /TDEPV/ Mq,MCl 

COMMON /MODES/ EV( 150 ,150 ) ,FREQ( 150 ) 

DO 10 1-1*3 
10 MC1(I)=0.0 

DO 40 1=1. N 
DO 45 J=I*3 
MQ( J , I )=MASS( I )*»Q( J , I ) 

45 MCK J)=MC1( J)«MQ( J,I) 

40 CONTINUE 

C COMPUTE "A" NEGLECTING TIME DEPENOENT( DEFORMATION DEPENDENT) TERMS 
C 

DO 50 J=l,3 
JP3=J*3 
DO 50 1=1,3 

50 A(I,JP3)=A12(I,J) 

DO 60 1=1*3 
IP3=I+3 
DO 60 J=l,3 
JP3=J^3 

60 A(IP3,JP3)= INERTl ( I , J ) 

DO 70 J=1*NT 

JP6=J«6 

DO 70 1=1*3 . 

IP3=I*3 

70 A(IP3,JP6)=A23(I,J) 

C 

C ADD IN FIRST ORDER DEFORMATION TERMS 
C 

CALL SKEU(MC1,UM) 

DO 80 J=l*3 
JP3=J*3 
DO 80 1=1*3 

80 A(I,JP3)=A(I*JP3)-WM(I.J) 

C FIRST ORDER DEFORMATION TERMS IN INERTIA MATRIX - "INERT2” 

C 

DO 85 1=1,3 
DO 85 J=l,3 
85 INERT2( I,J)=0.0 


00039300 

00039400 

00039500 

00039600 

00039700 

00039800 

00039900 

00040000 

00040100 

00040200 

00040300 

00040400 

00040500 

00040600 

00040700 

00040800 

00040900 

00041000 

00041100 

00041200 

00041300 

00041400 

00041500 

00041600 

00041700 

00041800 

00041900 

00042000 

00042100 

00042200 

00042300 

00042400 

00042500 

00042600 

00042700 

00042800 

00042900 

00043000 

00043100 

00043200 

00043300 

00043400 

00043500 

00043600 

00043700 

00043800 

00043900 

00044000 

00044100 

00044200 

00044300 

00044400 

00044500 
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ISN 0040 

ISN 004X 

ISN 0042 

ISN 0043 100 

ISN 0044 

ISN 0045 

ISN 0046 

ISN 0047 

ISN 0049 

ISN 0050 

ISN 0051 110 

ISN 0052 90 

ISN 0053 
ISN 0054 

ISN 0055 120 

ISN 0056 
ISN 0057 
ISN 0058 


GNNASS 


OS/360 FORTRAN H EXTENDED 


ISN 0059 

130 

ISN 0060 
ISN 0061 
ISN 0062 
ISN 0063 

140 

ISN 0064 
ISN 0065 
ISN 0066 
ISN 0067 
ISN 0069 
ISN 0070 

150 

ISN 0071 
ISN 0072 
ISN 0073 
ISN 0074 
ISN 0075 

160 

ISN 0076 
ISN 0077 
ISN 0078 

170 

ISN 0079 
ISN 0060 
ISN 0081 
ISN 0082 
ISN 0083 
ISN 0084 

165 

ISN 0085 

180 

ISN 0086 
ISN 0087 
ISN 0088 

C 

c / 

C 

190 


ISN 0089 
ISN 0090 
ISN 0091 
ISN 0092 
ISN 0093 
ISM 0094 
ISN 0095 
ISN 0097 


DO 90 L=l>3 
SUT1=0.0 
DO 100 I=1»N 
suii=suti+riRm Lfi) 

SUM=2.0«SU« 

LL=L*1 

DO 110 II’lt2 
IF(LL .GT, 3) LL=1 
INERT2( LL r LL)=INEHT2( LL » LL )4S«1 
LLSLL41 
CONTINUE 
CONTINUE 
SUM=0,0 
DO 120 I-1»N 

SUM=SUM- (tWt1(l.H<‘QC2.I),l*H(2.U»<«l.Il> 
INERT2(1.2)=Sm 
SUNsO.O 
DO 130 1=1. N 

SUri=SUM-( MRm 1 .1 3 .1)+riRM( 3.I)*»R( 1 . 1 ) ) 

INERT2(1.3)=SUtt 

sun=o,o 

DO 140 1=1 iN 

SUM=SUM-(MRM(2.I)»Q(3.I)^MRM( 3.I)»QC2.I) ) 
lNERT2t2.3)=SUM 
DO 150 1=1*3 
DO 150 J=li3 
IF(I ,LE. J) SO TO 150 
INERTEt I , J )=INERT2t J .1 ) 

CONTINUE 
DO 160 1=1.3 
IP3=I*3 
DO 160 J=l*3 
JP3=J*3 

A( 1P3 , JP3 )=AC IP3. JP3 )4IHERT2( I. J ) 

DO 170 1=1 iN 
L=3»I-2 

CALL SKEH(MQ(1.I).A23Q(1.L)) 

DO 180 1=1.3 
IP3=I+3 
DO 180 J=1*NT 
JP6=J«6 
DO 185 L=1.N3 

A(IP3* JP8 )=A( IP3» JP6)'fA23QC I.L)»EVC L. J) 
CONTINUE 

ADO IN SECOh© ORDER DEFORMATION TERMS - ''INEHT3'* 

DO 190 1=1.3 
DO 190 J=1.3 
INERT3(I>J)=0.0 
DO 200 L=l»3 
SUM=0.0 
DO 210 1=1. N 

210 SUt1=SUT1*MQ( L.I )*Qt L.I) 

LL=L*1 

DO 220 11=1.2 

IF(LL ,GT, 3) LL=1 

INERT3( LL,LL) = INERT3( LL.LD4SUM 


DATE 82.246/12.21.50 

00044600 

00044700 

00044800 

00044900 
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00045200 

00045300 
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00046000 

00046100 
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00046300 

00046400 

00046500 
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00046600 

00046900 
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00047100 

00047200 

00047300 

00047400 

00047500 

00047600 

00047700 

00047600 

00047900 

00046000 

00048100 

00048200 

00046300 

00046400 

00046500 

00048600 

00048700 

00046800 

00046900 

00049000 

00049100 

00049200 

00049300 

00049400 
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00049600 

00049700 
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00049900 

00050000 
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00050200 

00050300 
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DATE 82.2<»6/12.21.50 


PAGE 


ISN 0096 


LL=LL*1 


00050400 

ISN 0099 

220 

CONTINUE 


00050500 

ISN 0100 

200 

CONTINUE 


00050600 

ISN 0101 


SUM=0.0 


00050700 

ISN 0102 


DO 240 Isl,N 


00050000 

ISN' 0103 

240 

SUM=SUM-MQ( 1 , I )»Q( 2 , I ) 


00050900 

ISN 010^ 


INEHT3(1»2I=SUM 


00051000 

ISN 0105 


SUM=0 . 0 


00051100 

ISN 0106 


DO 250 I=:l>N 


00051200 

ISN 0107 

250 

SUM=SUM-MQC 1 , 1 )*Q( 3 » I ) 


00051300 

ISN 0106 


INERT3(1,3)=SUM 


00051400 

ISN 0109 


SUM=0.0 


00051500 

ISN 0110 


DO 260 I-1>N 


00051600 

ISN 0111 

260 

SUM=SUM-MQ ( 2 . I ) *Q( 3 , I ) 


00051700 

ISN 0112 


INERT3C2»3)=SUM 


00051800 

ISN 0113 


DO 270 I'l>3 


00051900 

ISN 0114 


DO 270 J-I.3 


00052000 

ISN 0115 


1F(I .LE. J) GO TO 270 


00052100 

ISN 0117 


INERT3( I , J )=INERT3( J , I) 


00052200 

ISN 0116 

270 

CONTINUE 


00052300 

ISN 0119 


DO 260 I=l>3 


00052400 

ISN 0120 


IP3=I43 


00052500 

ISN 0121 


DO 260 J=l»3 


00052600 

ISN 0122 


JP3=J43 


00052700 

ISH 0123 

260 

A(IP3,JP3)=A(IP3»JP3)+INERT3(I»J) 


00052600 


c 



00052900 

ISN 0124 

302 

FORMAT! 1H0,//»1H ,10X> ' A MATRIX* ) 


00053000 

ISN 0125 

301 

FORMAT(lHO,2X,lS(F7.2ilX)) 


00053100 

ISN 0126 


RETURN 


00053200 

ISH 0127 


END 


00053300 


*K)PTION3 IN EFFECT»NArtE(MAIN) NOOPTIMIZE LINECWJMTC 60 ) SI2E{03B^) AUTODBL(NONE ) 

^OPTIONS IN EFFECT>»SDURCE EBCDIC NOLIST NOOECK NOOBJECT NOttAP NOFORfUT GOSTMT NOXREF NOALC NOANSF TERH IBM FU6(I) 
•STATISTICS* SOURCE STATEMENTS = 126 » PROGRAM SIZE = 6156 > SUBPROGRAM NAME =6NMAS3 

•STATISTICS* NO DIAGNOSTICS GENERATED 


«»«««« end of COMPILATION •<hhhi* 


256K BYTES OF CORE NOT USED 


LEVEL Z.3.0 CJUHE 7«) M.2W./12.21.5E 

REQUESTED OPTlOHSt nOOBJ .TERM. .NOXREF . .H«UP. . .NAMEiriAIHI . AOCNONE ) ,OPT( 0 ) . . .FLAGd ) ,SIZE( 38<kK) ,LCC 60 > , 

.« .„.a. ■“ 


ISN 0002 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ISN 0003 


ISN OOOA 


ISN 0005 


ISN 0006 


ISN 0007 


ISN 0008 


ISN 0009 


ISN 0010 


ISN 0011 


ISN 0012 

10 

ISN 0013 


ISN OOIA 


ISN 0015 

20 

ISN 0016 


ISN 0017 

30 

ISN 0018 


ISN 0019 


ISN 0020 

50 

ISN 0021 


ISN 0022 


ISN 0023 

60 

ISN 00 2 A 

AO 

ISN 0025 


ISN 0026 


ISN 0027 

70 

ISN 0026 


ISN 0029 


ISN 0030 


ISN 0031 


ISN 0032 


ISN 0033 

85 

ISN 003A 

80 

ISN 0035 


ISN 0036 


ISN 0037 


ISM 0038 


ISN 0039 

95 


SUBROUTINE EXTF( FO ,TAU0 »FP ) 

TUTS SUBROUTINE ASSEMBLES THE FORCE VECTOR "F" IN THE 
JotiohIquatiqhs and is partitioned *s= ™ 

TRANSLATION, FORCES FOR BODY ROTATION, AHO FORCES FOR 
PARTICLE TRANSLATION. 

INPUT TO SUBROUTINE 

FO - EXTERNAL FORCE ON MO (AT 

TAUO - EXTERNAL TORQUE AT LOCATION OF BODY FRAME ORIGIN, 

FP - VECTOR OF EXTERNAL FORCES ON PARTICLES (L>2»..*iN. 

(ALL VECTORS EXPRESSED IN BODY FRAME) 

IMPLICIT PEAL»B( A-M*0“Z) 

OIMENSION^FO( 3)iTAU0(3)»FP(3»50)*Q(3*50)»RM( 3»S0) »MASS(50 

1 ^mSi.cmatJs^ 

2 OMEGA( 3 ) .R( 3 1 A THETA( 3 ) ,WV2( 150 ) *PHI( 150,150 ) ,WS( 150 ) 
COMMON /STATE/ R, THETA, Q,UVW, OMEGA. QOOT 

COfTON /CONST/ TMiMA5S,MCO,MRM,CMAT,M,N3,H3P6»HT,KTP6,NO 

COMMON /FORCE/ F 
COMMON /MODES/ PHI,W3 
DO 10 J>1,3 
F( J)-F0( J) 

00 20 1=1, N 
DO 20 J=l,3 
F( J)=F( J)*FP( J,I) 

DO 30 J=l,3 

SUM( J)=0.0 

DO AO 1=1 »N 

DO 50 J=l,3 

UV( J)=RM( J,I)«^Q( J»I) 

CALL CR0SS(WV,FP(1»I).WV1) 

DO 60 J=l,3 

SUM( J )=SUM( J )*HV1( J ) 

CONTINUE 
DO 70 J=l,3 
JP3=J*3 

FI JP3 )=TAU0( J )*SUM( J ) 

L=1 

DO 80 1=1 *N 
DO 65 J=l,3 
WV2( L)=FP( J.I) 

L=L+1 
CONTINUE 
CONTINUE 
00 90 1=1, NT 
IP6=I+6 
F(IP6 ) = 0.0 
DO 95 L=liN3 

F( IP6 )=F( IP6 )>PHI( L,I)»VIV2( L) 


0OO53AOO 
00053500 
00053600 
00053700 
00053800 
00053900 
0005A0OO 
0005A100 
0005A200 
0005A3DO 
0005AA0O 
0005A500 
0005A600 
0005A700 
0005A800 
0005A900 
),MC0( 3), 00055000 
),UVM( 3), 00055100 
00055200 
00055300 
00055AOO 
00055500 
00055600 
00055700 
00055800 
00055900 
00056000 
00056100 
00056200 
00056300 
00056AOO 
00056500 
00056600 
00056700 
00056800 
00056900 
00057000 
00057100 
00057200 
00057300 
00057A0O 
00057500 
00057600 
00057700 
00057800 
00057900 
00058000 
00058100 
00058200 
00058300 
00058A0Q 
00058500 
00056600 
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FLAG(I) 


52 


LEVEL 2.3,0 (JUNE 78) 


EXTF 
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DATE 82.246/12.21.52 


PAGE 


ISN 0040 90 CONTINUE 00058700 

ISN 0041 RETURN 00058800 

ISN 0042 END 00058900 

^OPTIONS IN EFFECT»NAME(I1AIN) NOOPTIMIZE LINECOUNT( 6 0 ) SIZE(0384K) AUTOOBL(NONE ) 

•OPTIONS IN EFFECT»SOURCE EBCDIC NOLIST NOOECK NOOBJECT NOMAP NOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM FLAG(I) 
•STATISTICS* SOURCE STATEMENTS = 41, PROGRAM SIZE = 2734, SUBPROGRAM NAME = EXTF 

•STATISTICS* NO DIAGNOSTICS GENERATED 


END OF COMPILATION ****** 


280K BYTES OF CORE NOT USED 


LEVEL 2.3.0 (JUNE 

requested OPTIONS! NOOBJ.TERMi .NOXREF » pNOTUP» 


OS/360 FORTPAN H EXTENDED DATE 62. 246/12. 2 X. 54 

, ,NAHE( MAIN ) » AD( NONE ) I OPT( 0 J , , » FLAG( I) f SIZE ( 364K ) , LC( 60 ) » 


OPTIONS IN effect: *«*“= HOAHSF TERN IBN 


ISN 0002 


ISN 0003 
ISN 0004 
ISN 0005 


ISN 0006 
ISN 0007 
ISN 0006 
ISN 0009 
ISM 0010 
ISN 0011 
ISN 0012 


c 

c 

c 


ISN 0013 
ISN 0014 
ISN 0015 
ISN 0016 
ISN 0017 

20 

ISN 0016 
ISN 0019 
ISN 0020 

30 

ISN 0021 
ISN 0022 
ISN 0023 

40 

ISN 0024 
ISN 0025 
ISN 0026 
ISN 0027 

42 

ISN 0026 
ISN 0029 
ISN 0030 

44 

ISN 0031 
ISN 0032 
ISN 0033 
ISM 0034 

55 

ISN 0035 
ISN 0036 
ISN 0037 
ISN 0036 

60 

ISN 0039 
ISN 0040 
ISN 0041 
ISN 0042 

70 

ISN 0043 

50 

ISN 0044 



SUBROUTINE NLKT 

TUTS subroutine CALCULATES THE HON-LINEAR KINEMATIC TERK^ IN THE 
^SnoSiSSJ^I^NS ^SE TERMS ARE ASSEMBLED INTO THE VECTOR -U». 

IMPLICIT REAL*6(A-H.0-Z) 

SlHE;^lS^l"M;fs? ^ *CMAT( 3. 3 NMQC 3.50 ) pMCK 3 ) , 

l”Son3?yoSf3l:ONEGA,I,.UT<3..UR<I..^^^^ 

2 MV2(3).WV3(3).WV4t3),WMll3p3).qi3t50)pR(3).THETA(3).M0Sl3). 

T Tor 1 .S( 3 ) »MS( 150 ) iPHK 150»150 ) 

C^ON /CONST/ TM.MASS.MCO.MRM, CHAT •N.N3.N3P6, NT.NTP6.no 

COMMON /TOEPV/ MQ.MCl 

COMMON /STATE/ R .THETA, q.UVW.OMESA.QOOT 
COMMON /FICFRC/ U 

COMMON /RIGID/ MOS.IO.S , 

COMMON /MODES/ PHI, MS 

EQUIVALENCEt U1 1 ) ,UT( 1 ) ) , < UJ 4 ) ,UR( 1 ) ) 

CALCUUTE DEFORMATION INDEPENDENT TERMS 
CALL CROSS(UVW. OMEGA, HVIJ * 

SUM=0ME6A( 1 )**MC0( 1 UOMEGAt 2 MH1C0( 2)^OMEGA( 3 )*»MC0C 3) 

OMSfOMEGAt 1 l**2+OM£GA( 2 )*Ht240MEGA( 3)*»2 
DO 20 J=l,3 

UT( J )sTM*UVl ( J )-SUM»OMESAI J )40MS*MCO( J ) 

CALL CfiO5S(MC0,WVl,WV2) 

MV3? J >=CMAT( J , 1 )»K)MEGA( 1 )+CMAT( J ,2 )IK3MEGA( 2 )>CMAT( J » 3 )*OMEGA( 3 ) 
CALL CROSS ( OMEGA. MV3,HV4) 

DO 40 J=1.3 

^=WWA( i jis): I )*01EGA( 2 )»S( 2 )40MEBA< 3 )»S( 3 1 
CALL CROSS* MO S, OMEGA, HV2) 

HV3( J )=IO( J ,1 )»OMEGA* 1 )*I0( J,2 )«K5ME6A( 2 )*I0C J,3 l*OMEGA( 3 ) 

CALL CROSS* ONEGA, MV3.WV4) 

DO 44 J>I,3 

UR( J )=UR< J )-SUM»WV2U )-HV4C J ) 

DO 50 1=1. N 
L=3*I-2 
DO 55 J=li3 

i‘ n*«1E6A( , )««M( 2.1 

00 60 J=l,3 
WV3( J 1 =SUtt»OMEGA( J ) 

WV4< J )=OMS»MRM( J , I) 

LL=L 

DO 70 J=1.3 

UV( LLJ=MV2( J )tWV 3( J UWV4t J I 
LL=LL+1 
CONTINUE 
DO 74 1=1. NT 
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00062100 
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LEVEL 2.3.0 (JUNE 76) 


NLKT 
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PAGE 2 


ISN 

0045 


IP6=I^6 

00064300 

ISN 

0046 


U(IP6)=0.0 

00064400 

ISH 

0047 


DO 74 L=1,N3 

00064500 

ISN 

0048 

74 

U(IP6)=U(IP6)^PHI( L,I)»UV( L) 

00064600 



c 


00064700 



c 

CALCULATE FIRST ORDER DEFORMATION DEPENDENT TERMS 

00064800 



c 


00064900 

ISN 

0049 


DO 60 J=l.3 

00065000 

ISN 

0050 

60 

HV2( J]=0.0 

00065100 

ISN 

0051 


DO 90 1=1 »N 

00065200 

ISN 

0C52 


DO 95 J=l>3 

00065300 

ISN 

0053 

95 

WV2( J )=WV2( J )+MASS( I )*QDOTC J ,I ) 

00065400 

ISN 

0054 

90 

CONTINUE 

00065500 

ISN 

0055 


CALL CROSS( OMEGA, WV2,WV3) 

00065600 

ISN 

0056 


SUM=OMEGA( 1 )*»MClt 1) ♦OMEGA ( 2 )M1C1C 2 )*OMEGA( 3 )**«C1( 3 ) 

00065700 

ISN 

0057 


DO 100 J=l,3 

00065800 

ISN 

0058 

100 

UT( J )=UT( J )-2 . 0**MV3( J )-SUM»»OMEGA( J )40MS»MC1( J ) 

00065900 

ISN 

0059 


CALL CROSS(MCl,WVl,WV2) 

00066000 

ISN 

0060 


DO 110 1=1,3 

00066100 

ISN 

0061 


DO 110 J=lt3 

00066200 

ISN 

0062 

110 

UM1(I,J)=0.0 

00066300 

ISN 

0063 


00 120 1=1 »N 

00066400 

ISN 

0064 


DO 125 J=I*3 

00066500 

ISN 

0065 

125 

HMK J,J)=MM1( J,J)+2.0«MRM(J,I)»Q(J,I) 

00066600 

ISH 

0066 


W11Cl.2)=NMl(l,2)+MRM(l,I)*»Q(2,I)*MRM(2,I)*Q(l.I) 

00066700 

ISN 

0067 


WMl(1.3)=WMl(l,3)+MRM(l,I)*Q(3,I)*MRM(3,IH«J(liI) 

00066800 

ISH 

0068 


HMl(2,3)=WMlt 2,3)*MRM(2,I)*Q( 3,I)4MRM(3,I)»Q(2,I) 

00066900 

ISN 

0069 

120 

COtfTINUE 

00067000 

ISN 

0070 


DO 130 1=1,3 

00067100 

ISN 

0071 


DO 130 J=l,3 

00067200 

ISN 

0072 


IF(I .LE. J) 60 TO 130 

00067300 

ISN 

0074 


NM1(I,J)=MM1( J,ll 

00067400 

ISN 

0075 

130 

CONTINUE 

00067500 

ISN 

0076 


DO 140 J=l,3 

00067600 

ISN 

0077 

140 

HV3( J )=WM1( J , 1 )*OMESA( 1 J4WM1( J ,2 )»OMEGA( 2 )*WM1( J , 3 )*0ME6A( 3 ) 

00067700 

ISN 

0078 


CALL CROSS! OMEGA, WV3,WV4) 

00067800 

ISH 

0079 


SUM=0.0 

00067900 

ISN 

0060 


DO 150 1=1, N 

00066000 

ISN 

0081 

150 

SUM=SUM4MRM( 1 , I )«QDOTt 1 , I )4MRM( 2 ,I )«OOOT( 2,1) ♦MRM( 3 , I )*Q0OT( 3,1) 

00068100 

ISN 

0082 


DO 160 J=li3 

00068200 

ISN 

0083 

160 

UVl ( J )=SUM»0 MEGA( j ) 

00068300 

ISN 

0084 


DO 170 J=l,3 

00068400 

ISN 

0085 

170 

HV3( J)=0.0 

00068500 

ISN 

0086 


DO 180 1=1, N 

00068600 

ISN 

0087 


SUM=MRM( 1 , r )*OMESAt 1 )*MRM( 2 , 1)»0MEGA( 2 )*«?M( 3 , 1 )*0ME6A( 3 ) 

00068700 

ISN 

0088 


DO 190 J=l,3 

00066800 

ISN 

0089 

190 

WV3< J )=WV3( J )4SUM»*Q00T( J, I ) 

00068900 

ISN 

0090 

180 

CONTINUE 

00069000 

ISN 

0091 


DO 200 J=l,3 

00069100 

ISN 

0092 

200 

UR( J)=UR(J)+HV2tJ)+WV4( J)^2.»(HV3(J)-WV1( J) ) 

00069200 

ISN 

0093 


DO 300 1=1, N 

00069300 

ISH 

0094 


L=3*I-2 

00069400 

ISN 

0095 


CALL CROSS ( OMEGA, QOOTCl, I ),WV1) 

00069500 

ISN 

0096 


DO 310 J=l,3 

00069600 

ISN 

0097 

310 

WV2( J )=MASS( I )»WV1< J ) 

00069700 

ISN 

0098 


StR1=OMEGA( 1 )«MQ( 1 ,1 )^OMEGA( 2 )**MQ( 2,1 )+OMEGA( 3)«MQ( 3,1 ) 

00069800 

ISN 

0099 


DO 320 J=l,3 

00069900 

ISN 

0100 


WV3( J ) =SUM»CMEGA( J ) 

00070000 


55 


LEVEL 2.3.0 (JUNE 7S) 


NLKT 


OS/3AO FORTRAN H EXTENDED 


ISN 0101 
ISN 0102 
ISN 0103 
ISN 0104 
ISN 0105 
ISN 0106 
ISN 0107 
ISN 0108 
ISN 0109 
ISN 0110 
ISN 0111 


ISN 0112 
ISH 0113 
ISN 0114 
ISN 0115 
ISN 0116 
ISN 0117 
ISN 0118 
ISN 0119 
ISN 0120 
ISN 0121 
ISN 0122 
ISN 0123 
ISN 0124 
ISN 0126 
ISN 0127 
ISN 0128 
ISN 0129 
ISN 0130 
ISN 0131 
ISN 0132 
ISN 0133 
ISN 0134 
ISN 0135 
ISN 0136 
ISN 0137 
ISN 0138 
ISN 0139 
ISN 0140 
ISN 0141 
ISN 0142 
ISN 0143 
ISN 0144 
ISN 0145 


320 WV4lJ>=0MS»nq( J»I) 

LL=L 

Sv( LU =UVl i.L ) -2 . »HV2 ( J ) -HV3( J J ) 

330 LL=LU*1 
300 CONTINUE 

00 307 1=1 »NT 
IP6=I+6 
UCIP6)=0.0 
00 307 L=liN3 

307 U( IP6 J=U( IP6 )^PHIJ L^D^UVl L) 

C CALCULATE SECOND ORDER DEFORHATION DEPENDENT TERMS 


311 


331 


321 


340 

350 


360 


370 


385 

380 

390 


DO 311 1=1.3 
DO 311 J=1.3 
wm(i.J)=o.o 
DO 321 1=1. N 

DO 331 J=l»3 . 

UMK J. J )=WM1( J. J.I J* J J 

Mril(1.2)=WMltl»2»+MQtl.I)»<3t2,I) 

wii(i.3)=w«ui.3uriQti.i)*«i3.n 

WMK 2.3)=WMlt 2»3)*MQl 2.1 3*11 

CONTINUE 

DO 340 1=1.3 

DO 340 J=l*3 

IF(I .LE. J) 60 TO 340 

wMia.J)=MMi(j,n 

CONTINUE 

£5i?j%wni( 1»aiEGAt 1 )*W11( J.2 2 J.3»*0MESAI 31 

CALL CROSSIOMEGA.WVI.HVE) 

sun=o.o 

nSa . I ).<W0Ta . I . 2 . 1 )«W0T( 2 .1 
DO 370 J=1.3 
WV3t J 1 =SUM»0ME6A( J ) 

UV4U)=0.0 

5u«=OMEGi( n-rwi 1 ,1 ) ♦0«EGA( 2 )1im 2.1 1^E6A< 31fMI 3,11 
DO 385 J=1.3 

mv4( j )=wv4( j )4sun»<aooT{ j. I ) 

CONTINUE 

RETURN 
END 


DATE 82.246/12.21.54 


00070100 

00070200 

00070300 

00070400 

00070500 

00070600 

00070700 

00070800 

00070900 

00071000 

CC071100 

00071200 

00071300 

00071400 

00071500 

00071600 

00071700 

00071800 

00071900 

00072000 

00072100 

00072200 

00072300 

00072400 

00072500 

00072600 

00072700 

00072800 

00072900 

00073000 

00073100 

00073200 

00073300 

00073400 

00073500 

00073600 

00073700 

00073800 

00073900 

00074000 

00074100 

00074200 

00074300 

00074400 

00074500 

00074600 

00074700 


PAGE 


ISN — 

.OPTIONS IN EFrECT.NAnE<«AIN. NOOPTINIZE LINEC<XKT(60 . SIZEC.3SAK. AUTOOBU NONE I 

.OPTIO,« IN EFFECT.SOURCE EBCOIC HOLIST NOOECK NOOBJECT NONAP NOFORHAT GOSTNT NOXBEF NOALC NO^F TEPN 


..STATISTICS* SOURCE STATEMENTS » 

:»STATISTICS* NO DIAGNOSTICS GENERATED 
»**.!** END OF COMPILATION ****** 
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244K BTTES OF CORE NOT USED 


56 



LEVEL 2.3.0 (JUNE 78) OS/360 FORTRAN H EXTENDED DATE 62.2W12.21.56 

REQUESTED OPTIONS: NOOBJ »TERM» tNOXREF , »NOMAP, , rNAME(rtAIN)fAD(NONE )tOPT( 0)i t >FLA6( I ) »SIZE( 38^K ) » tC( 60 ) » 

OPTIONS IN effect: NAME (MAIN) N00PTIMI2E LINECOUNTC 60 ) SIZE(038^K) AUTOOBU NONE ) 

SOURCE EBCDIC NOLIST NODECK NOOBJECTT NOMAP NOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM 


ISN 

0002 


SUBROUTINE SOLVE 




00076800 



c 





00076900 



c 

THIS SUBROUTINE ASSEMBLES THE EQUATIONS OF MOTION IN 

FIRST 

ORDER 

00075000 



c 

FORM AND SOLVES THE SET OF SIMULTANEOUS DIFFERENTIAL 

EQUATIONS 

00075100 



c 

(SIZE 6«N^12) 




00075200 

ISN 

0003 


IMPLICIT REAL»8(A-H,0-Z) 




00075300 

ISN 

000^ 


REAL«8 MASS,MC0,MRM,INERT1 




00075600 

ISN 

0005 


DIMENSION MASS150 ) ,MC0( 3 ) »MRM( 3i50 ) tCMAT( 3*3 ) 

,A(156 

,156). 

A12(3, 

3), 00075500 




1 INERT1(3.3)*A23(3.150),Q(3.50),QDOT(3,50),UVW(3)*OMEGA(3). 

00075600 




2 Rt3)*THETA(3)*F(156),U( 156),R16(3*3),GAMA(3, 

3), 



00075700 




3 B( 156 ) , YDOT( 312 ) * Y( 312 ) ,WV1( 3 ) ,HV2( 3) *WV3( 156 ) . 



00075800 




6 QV(1501,QOOTV( 150 )*AV( 12266) 




00075900 

ISN 

0006 


EQUIVALENCE ( QV( 1 ) .Q( 1 » 1 ) ) . ( QOOTV( 1 ) ,QOOT( 1,1)) 



00076000 

ISN 

0007 


COMMON /CONST/ TM, MASS, MC 0 ,MRM, CHAT, N,N3,N3P6 

,NT,NTP6,NO 


00076100 

ISN 

0008 


COMMON /AMAT/ A,A12 , INERTl , A23 




00076200 

ISN 

0009 


COMMON /STATE/ R , THETA, Q,UVM, OMEGA, (JOOT 




00076300 

ISN 

0010 


COMMON /FORCE/ F 




00076600 

ISN 

0011 


COMMON /FICFRC/ U 




00076500 

ISN 

0012 


COMMON /TIME/ OT,TSTOP,DTP,DTG 




00076600 

ISN 

0013 


COMMON /MODES/ PHK 150 , 150 ) ,WS( 150 ) 




00076700 

ISN 

0016 


COMMON /MODCO/ ETAU50 ) ,ETAO( 150 ) 




00076600 

ISN 

0015 


DATA GAMA/e^O.O.l.O/rlPASS/O/ 




00076900 



c 





00077000 



c 

CALCULATE R16 -TRANSFORMATION FROM BOOT FRAME TO INERTIAL FRAME 

00077100 



c 





00077200 

ISN 

0016 


Sl=OSIN(THETA(D) 




00077300 

ISN 

0017 


C1=OCOS(THETA(D) 




00077600 

ISN 

0018 


S2=OSIN(THETA(2)) 




00077500 

ISN 

0019 


C2=0C0S(THETA(2)) 




00077600 

ISN 

0020 


S3=DSIN(THETA(3) ) 




00077700 

ISN 

0021 


C3=0C0S(THETA(3)) 




00077800 

ISN 

0022 


R16(1,1)=C2»C3 




00077900 

ISN 

0023 


R16C1,2)=-C2«S3 




00078000 

ISN 

0026 


R16(1,3)=S2 




00078100 

ISN 

0025 


R16( 2 , 1 )=Cl»S3*Sl*tS2*C3 




00076200 

ISN 

0026 


R16(2,2)=Cl»C3-Sl»tS2»S3 




00078300 

ISN 

0027 


R16(2,3)=-S1»C2 




00078600 

ISN 

0028 


R16( 3,1 )=S1»S3-C1*S2*C3 




00076500 

ISN 

0029 


R16( 3,2 )=S1*»C3*CI*S2«S3 




00078600 

ISN 

0030 


R16(3,3)=C1*»C2 




00078700 



c 





00078800 



c 

CALCUUTE **GAMA“ - TRANSFORMS ANGULAR VELOCITY 

TO ATTITUDE 

RATES 

00078900 



c 





00079000 

ISN 

0031 


GAMA(1,1)=C3/C2 




00079100 

ISN 

0032 


GAMA(1,2)=-S3/C2 




00079200 

ISN 

0033 


GAMA(2,1 )=S3 




00079300 

ISN 

0036 


GAMA(2,2)=C3 




00079600 

ISN 

0035 


T2=S2/C2 




00079500 

ISN 

0036 


GAMAC3,1)=-C3*T2 




00079600 

ISN 

0037 


GAMA ( 3,2 )=S3»T2 




00079700 

ISN 

0038 


DO 30 1=1,6 




00079600 

ISM 

0039 

30 

B(I)=F( I)*U( I) 




00079900 

IS)( 

0060 


DO 60 1=1, NT 




00080000 
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SOLVE 


OS/360 FORTRAN H EXTENDED 


DATE 62,266/12.21.56 


PAGE 2 


LEVEL 2.3.0 (JUNE 78) 


ISN 0061 
ISN 0062 


ISN 0062 
ISN 0063 
ISN 0066 
ISN 0065 
ISN 0066 
ISN 0067 
ISN 0066 
ISN 0069 
ISN 0070 
ISN 0071 


60 

C 

c 

c 


ISN 

0063 


ISN 

0066 


ISN 

0065 


ISN 

0066 


ISN 

0067 


ISN 

0066 

300 

ISN 

0069 


ISN 

0050 


ISN 

0052 


ISN 

0053 


ISN 

0056 

61 

ISN 

0055 

60 

ISN 

0057 

C 



c 



c 

ISN 

0056 


ISN 

0059 


ISN 

0060 


ISN 

0061 

70 


60 


90 


100 

c 

c 

c 


ISN 

0072 

105 

ISN 

0073 


ISN 

0076 

110 

ISN 

0075 


ISN 

0076 


ISN 

0077 

120 

ISN 

0076 


ISN 

0079 

130 

ISN 

0060 


ISN 

0081 

160 



C 



c 



c 

ISN 

0062 


ISN 

0063 


ISN 

0086 


ISN 

0065 

150 

ISN 

0036 



ISN 0067 


IP6sl^6 

B( IP6 )=F( IP6 )4U( IP6 )-MS( I )*ETA( I ) 

STORE UPPER TRIANGLE OF “A’* IN *'AV** 

L=l 

00 300 J=l,NTP6 
00 300 I^lfJ 
AV(L)sA(IiJ) 

L«L*1 

SlL^LEQT 1P( AVil»NTP6.B.156»0*01*D2»IEH ) 

IFdER .EQ. 0) 60 TO 60 
URITE(6»61) lER 

mSmTUHO.SX. -ERROR DETECTED BY IHSL tIBRARY ROUTINE "LEQTIP" 
IERROR C00E=*»I3) 

IFdPASS ,EQ. 1) GO TO 105 
IPASS=1 

SET INITIAL VALUE OF 

• DO 70 1=1*3 
IP3=I^3 
Yd)=Rd) 

YdP3)=TNETAd) 

DO 60 1=1 pNT 
Yd+6)=ETAd» 

DO 90 1=1*3 
LsHTP6^I 
LL=L*3 
Y( L)=UVM(I) 

Y(LL)=OrtEGAd> 

DO 100 1=1 iNT 

LS12+NT+I 

Y(L)=ETAOd) 

SET UP ’•YOOT" 

DO 120 1=1*3 
YDOT( I )=WV1( I ) 

YD0Td43)=WV2d) 

00 130 1=1 *NT 
YDOT(64l)=ETAD(I) 

DO 160 I=1.NTP6 
YOOT(MTP6 + I)=Bd) 

UPDATE VARIABLES IN STATE VECTOR 

CALL ODESLV(NO*Y*YDOT.DT) 

DO 150 1=1*3 
Rd) = Yd) 

THETAd)=Yd43) 

DO 160 1=1. NT 
FTAdl=Y<6•^I) 


00080100 

00060200 

00080300 

00060600 

00060500 

00060600 

00060700 

00080600 

00060900 

00061000 

00061100 

00081200 

00061300 

00061600 

00061500 

00061600 

00061700 

00081800 

00061900 

00062000 

00062100 

00062200 

00062300 

00062600 

00062500. 

00082600 

00082700 

00062600 

00062900 

00083000 

00063100 

00063200 

00083300 

00063600 

00063500 

00083600 

00063700 

00083600 

00063900 

00066000 

00066100 

00086200 

00066300 

00066600 

00066500 

00086600 

00066700 

00086600 

00066900 

00065000 

00065100 

00085200 

00065300 

00085600 

00085500 

00085600 

00065700 

00065600 


t 

c 

I 

n 

I 
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LEVEL 2,3.0 (JUNE 76) 


SOLVE 


OS/360 FORTRAN H EXTENDED 


DATE 82.2<*6/12.21.56 


PAGE 3 


ISN 

0088 


DO 170 I<1*3 

00085900 

ISN 

0089 


L=NTP6+I 

00086000 

ISN 

0090 


LL=L*3 

00086100 

ISN 

0091 


UVW(I)=Y(L) 

00086200 

ISN 

0092 

170 

OMEGA(I) = Y(LU 

00086300 

ISN 

0093 


DO 180 1=1. NT 

00086400 

ISN 

0094 

180 

ETA0(I)=YCNT>124I) 

00086500 


ISN 0095 
XSN 0096 
ISN 0097 
ISN 0098 
ISN 0099 
ISN 0100 
ISN 0101 
ISN 0102 
ISN 0103 


220 

200 


COMPUTE NEW VALUES FOR "Q” AND •'ROOT*' 

00 200 I=1>N3 
QV(I)=0.0 
QDOTV(IJ=0.0 
00 220 L-liNT 

<JVCI)=QV(I) + PHI(I,L)»ETA(L) 

qOOTV( I )=QOOTV( I )^PHI( I » L J»ETAO( L) 

CONTINUE 

RETURN 

END 


00086600 

00086700 

00086800 

00086900 

00087000 

00087100 

00087200 

00087300 

00087400 

00087500 

00087600 

00067700 


•OPTIONS IN EFFECT«NAf1E(MAIN) NOOPTIMIZE UNECOUMT(60) SIZE(0384K) AUTOOBUNONE) 


•OPTIONS IN EFFECT»SOURCE EBCDIC NOLIST NOOECK N008JECT NOMAP NOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM FLAG(I) 
•STATISTICS^ SOURCE STATEMENTS = 102 » PROGRAM SIZE s 107808 » SUBPROGRAM NAME = SOLVE 


•STATISTICS* NO DIAGNOSTICS GENERATED 

END OF COMPILATION <hhhhhi 256K BYTES OF CORE NOT USED 
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level E.3.. .JUNE 7») 

REQUESTED OPTIONS: NOOOJ ,TERH. .NOXREF. .NOHAP. . .NAHE(IIAIH) .AOCHOHE . .OPTC 0 . , . .FLAQ. I . .SIZEl 384K . . LCl 60 ) . 

OPTIONS IN effect: «,*LC HOAHSF TERN IBN 


I5N 0002 
ISN 0003 
ISN 0004 
ISN 0005 


ISN 0006 
ISN 0008 
ISN 0009 
ISN 0010 

ISN 0011 


c 

c 

c 

c 

c 

c 

c 

c 

c 


12 

10 

c 

c 

c 


SUBROUTINE OOESLV(NtYiDERIV»H ) 

skl'oerwInk;^^^ 

data IHTF/l/ . CX/O . 0/ ,C2/0 . 0/ . C3/0 . / 

•mis SUBROUTINE INTEGRATES TME FIRST ORDER STSTEN OF OTMHARY 

m”er^?^ ewations ■■dy/dt=deriv by tne aoans nethod 

USING THIRD ORDER DIFFERENCES. 

t VEC™'or?NniAt VALUES OH INPUT. "V IS OVERV.ITTEM 
WITH THE NEW SOLUTION 
H- STEP SIEE 

IF(N .LE. 312) GO TO 10 
WRITE(6»12) N 

fSr^KIHO.SX. -ERROR IN 

ISIZE exceeds DIMENSION SIZE OF ARRAYS ) 

GO TO( 1000* 2000 i3000 »4000 )»INTF 


FIRST CALL TO ROUTINE - EULER INTEGRATION 


ISN 

0012 

1000 

ISN 

0013 

20 

ISN 

0014 


ISN 

0015 

c 



C 1 



C 

ISN 

0016 

2000 

ISN 

0017 


ISN 

0018 


ISN 

0019 

30 

ISN 

0020 


ISN 

0021 


ISN 

0022 



DO 20 1=1 »N 

DERIVOa)=DERIVm 

INTF=2 

GO TO 5000 


SECOND CALL TO ROUTINE - FIRST ORDER DIFFERENCES 


DO 30 1=1 »N 

BDK 1 ,1 )=OERIV( I J-DERIVOl I ) 
B01tl»2)=BDllIil) 
DERIVO(I)=OERIV( I) 

Cl=.5 
INTF=3 
GO TO 5000 


ISN 0023 
ISM 0024 
ISN 0025 
ISN 0026 
ISN 0027 
ISN 0028 
ISN 0029 
ISN 0030 
ISN 0031 


ISN 0032 
ISN 0033 


C 

c 
c 

3000 


40 


c 
c 
c 

4000 


THIRD CALL TO ROUTINE - SECOND ORDER DIFFERENCES 


DO 40 l=liN 

BOU I»2 )=OERIV( I )-OEHlVO( I) 
BD2tI»l)=B01( 1 »2 )-B01( I»l) 
B02tI>2)=B02(Iil) 
DERIVO(I1=OERIV(I) 
BD1(I»1)=B01CI»2) 

INTF=4 
C2=5. 0/12.0 
60 TO 5000 


ADAMS METMM) HITH 3RD ORDER DIFFERENCES. 
DO 50 1=1>N 

BDK I »2 )=DERIV( 1 )-DERIVO( I) 


000B7B00 
00087900 
)00068000 
00088100 
00088200 
00088300 
00068400 
00088500 
00088600 
00088700 
00088800 
00088900 
00089000 
00089100 
00089200 
00069300 
00089400 
00089500 
00089600 
00069700 
00089600 
00089900 
00090000 
00090100 
00090200 
00090300 
00090400 
00090500 
00090600 
00090700 
00090800 
00090900 
00091000 
00091100 
00091200 
00091300 
00091400 
00091500 
00091600 
00091700 
00091800 
00091900 
00092000 
00092100 
00092200 
00092300 
00092400 
00092500 
00092600 
00092700 
00092800 
00092900 
00093000 


PAGE 1 


FLAG(I) 


60 



LEVEL 2.3.0 (JUNE 78) 


OOESLV 


03/360 FORTRAN H EXTENDED 


DATE 82. 246/12. 21. 58 


PAGE 


ISN 

0034 


B02(I,2)=BD1(I,2)-B01(I,1) 

0OO931OO 

ISN 

0035 


BD3(I)=&02(I,2)>6D2(I,1) 

00093200 

ISN 

0036 


DERIVO(I):DERIV(I) 

00093300 

ISN 

0037 


BD1CI,1)=BDI(I,2) 

00093400 

ISN 

0038 

50 

B02(I,1)=BD2(I,2) 

00093500 

ISN 

0039 


C3=3. 0/8.0 

00093600 

ISN 

0040 


60 TO 5000 

00093700 



C 


00093800 



c 

UPDATE VECTOR *Y' 

00093900 



C 


00094000 

ISN 

0041 

5000 

DO 60 Isl,N 

00094100 

ISN 

0042 

60 

Y( I)=Y( I ]«K*(DERIV( I)4^C1*BD1(I,2)«C2»B02(I.2)4C3*BD3( I) ) 

00094200 

ISN 

0043 


RETURN 

00094300 

ISN 

0044 


END 

00094400 

♦OPTIONS IN 

EFFECT*NAME(MAIN) NOOPTIMIZE LINECOUNTC60) SIZE(0384K) AUTODBL( NONE ) 



♦OPTIONS IN EFFECT*SOURCE EBCDIC NOLIST NCDECK NOOBJECT NOMAP HOFORNAT GOSTMT NOXREF NOALC NOANSF TERM IBM FLAG(I) 
♦STATISTICS* SOURCE STATEMENTS = 43 » PROGRAM SIZE = 16838, SUBPROGRAM NAME =OOESLV 

♦STATISTICS* NO DIAGNOSTICS GENERATED 

END OF COMPILATION *♦«*** 276K BYTES OF CORE NOT USED 


LEVEL 2.3,0 I JUNE 
REQUESTED OPTIONS: 


03/360 FORTRAN H EXTENDED 

NOOBJiTERHi »NOXREF» »NOI1APi , ,NAMEIHAIN) »AO(NONE ) »0PT( 0 ) 


DATE 82.2A6/12.21.59 
, , ,FU6l I ) , SIZE( 38<VK ) , LCl 60 ) » 


c.<«. .H .„.=T. “*»' ™" ■“ 


I5N 0002 
ISN 0003 
ISM 00 0<* 
ISH OOOS 
ISN 0006 


ISH 0007 
ISN 0008 
ISN 0009 
ISH 0010 
ISN 0011 
ISN 0012 
ISN 0013 
ISN 001^ 
ISN 0015 
ISN 0016 
ISN 0017 
ISN 0018 
ISN 0019 
ISN 0020 
ISN 0021 
ISN 0022 
ISN 0023 


ISN 0024 
ISN 0025 
ISN 0026 

ISN 0027 
ISH 0028 
ISN 0029 
ISH 0030 


J 


ISN 0032 
ISM 0033 
ISN 0035 
ISN 0036 
ISH 0037 
ISN 0038 
ISN 0040 


ISN 0042 
ISN 0043 
ISN 0044 


201 


10 


20 

30 

100 

110 


120 

130 

200 

202 


C 

c 

c 


c 

c 

c 


310 

C 

C 

c 


SUBROUTINE PRINT(T.F0 »TAU0,FP) 

IMPLICIT REAL»8(A-H»0-Z) 

REAL«8 MASS,MCO»MRM 

dimension' R( 3')?THETA( 3 ) ,QI 3.50 ) .UVM( j®5o’, 

1 MC0(3).MRM(3.50),C«ATl3,3).IUSS«5O»,F0t3».TAUO(3>,FPt3,SO), 

2 IPLOTI42)#PLTDAH100,20 J 

COMMON /STATE/R I THETA »Q>UVW, OMEGA fQDOT 

COMMON /CONST/ TM,MASS,MCO ,MRM,CMAT ,NfN3»N3P6 ,NT,NTP6 »NO 

COMMON /PLOTT/ PLTDAT»IPLOT»NP 

MRIT£(6>100) T 

WRITE(6f200) FOfTAUO 

DO 201 I=lfN 

NRITE(6»202) I.IFPC J*D > 

DO 10 1=1*3 

TD( I )=57. 29578**THETA( I ) 

00( I )=57.29578*OMEGA( I ) 

URITE(6*110) R.TD,UVW,00 
DO 20 I=1*N 

URITE(6il20) I*tQ(JiI>*J-I#3) 

DO 30 I=liN 

WRITE(6*130) I»(QOOTt JiI)*J=l*3) 

FORMAT(1HO*///*5X,'TIME=**F10.3,’ SEC|) 

FORMATtlH ,6 X*‘R=* f3C2X»1PE11.4), »/,4X* , 

1 3(2X.1PE11.4)*‘ DEG’,/.6X,*UVM=:*,3(2X*1PE11.4)»* FT/SEC */. 

2 4X, ‘OMEGAS* »3(2X»1PE11. 4)** 0E6/SEC*,/J 
FORMATtlH ,5X,‘Q**I3,‘=**3(2X.1PE11,4).* FT') 


FORMATtlH i2X, 'QDOT* ,I3**=* ,3(2X»1PE11.4), 

FORMATtlHO*aX,*FO=' .3(2X*1PE11.4)# LB */* 

1 6X**TAU0=' .312X,1PE11.4)*' FT LB*) 

FORMATtlH ,4X» *F' *13* '=* »3(2X*1PE11.4)* LB ) 
RETURN 

ENTRY GRAF(T) 

IFtlPLOT(l) .EQ. 0) GO TO 203 
STORE VARIABLES FOR PLOTTING 


FT/SEC* ) 


NP=NP*1 

IFtNP .6T. 100) 60 TO 203 
PLTOATtNP,l)=T 
DO 300 1=1*42 
NYslPLOT(I) 

IFtNY .EQ. 0) SO TO 203 
IFtNY .GT. 3) GO TO 310 

STORE INERTIAL POSITION 

PLTOATtNP,I*l)=R(NY) 

60 TO 300 

IFtNY .GT. 6) GO TO 320 
STORE ATTITUDE ANGLES IN DEGREES 


00094500 
00094600 
00094700 
00094800 
00094900 
00095000 
00095100 
00095200 
00095300 
00095400 
00095500 
00095600 
00095700 
00095800 
00095900 
00096000 
00096100 
00096200 
00096300 
00096400 
00096500 
00096600 
00096700 
00096800 
00096900 
00097000 
00097100 
00097200 
00097300 
00097400 
00097500 
00097600 
00097700 
00097800 
00097900 
00098000 
00098100 
00098200 
00098300 
00096400 
00098500 
00098600 
00098700 
00098800 
00098900 
■ 00099000 

00099100 
00099200 
00099300 
00099400 
00099500 
00099600 
00099700 
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FUGtl) 


62 



LEVEL 2.3.0 (JUNE 76) 


PRINT 


03/360 FORTRAN H EXTENOEO 


DATE 82.266/12.21.59 


PAGE 2 


ISN 

0046 


PLTD AT ( NP , I ♦ 1 ) =57 . 29578*THETA ( NY-: 

ISN 

0047 


60 TO 300 

ISN 

0048 

320 

C 

C 

IF (NY .GT. N3P6) GO TO 330 
STORE DEFORMATION COORDINATES 

ISN 

0050 

c 

L=NY-6 

ISN 

0051 


Il=l*L/3 

ISN 

0052 


I2=L-3*(I1-1) 

ISN 

0053 


IF(I2 .NE. 0) GO TO 321 

ISN 

0055 


11=11-1 

ISN 

0056 


12=3 

ISN 

0057 

321 

PLTDAT(NP,I+1)=Q( 12,11) 

ISN 

0058 


60 TO 300 

ISN 

0059 

330 

IF(NY .GT. (N3P64^3)) GO TO 340 


ISN 0061 
ISN 0062 
ISN 0063 


ISN 0065 
ISN 0066 


340 

C 

c 

c 


STORE TRANSLATIONAL VELOCITY 

PLTO AT( NP , H^l ) =UVW( NY-N3P6 ) 

GO TO 300 

IF(NY .GT. (N3P6^6)) GO TO 350 

STORE ANGULAR VELOCITY IN DEG. /SEC. 

PLTOAT( NP ,1+1 ) =57. 29578*0NE6A( NY-( N3*9 ) ) 
60 TO 300 

STORE DEFORMATION RATES 


00099800 

00099900 

00100000 

00100100 

00100200 

00100300 

00100400 

00100500 

00100600 

00100700 

00100800 

00100900 

00101000 

00101100 

00101200 

00101300 

00101400 

00101500 

00101600 

00101700 

00101800 

00101900 

00102000 

00102100 

00102200 

00102300 

00102400 

00102500 

00102600 


ISN 

0067 

350 

L=NY-(N3+12) 

00102700 

ISN 

0068 


11=14L/3 

00102800 

ISN 

0069 


I2=L-3*(I1-1) 

00102900 

ISN 

0070 


IF(I2 .NE. 0) 60 TO 351 

00103000 

ISN 

0072 


11=11-1 

00103100 

ISN 

0073 


12=3 

00103200 

ISN 

0074 

351 

PLTOAT(NP,I-H)=QOOT( 12,11) 

00103300 

ISN 

0075 

300 

CONTINUE 

00103400 

ISN 

0076 

203 

RETURN 

00103500 

ISN 

0077 


END 

00103600 


•OPTIONS IN EFFECT*NAME(MAIN) HOOPTIMIZE LINECOUNT( 60 ) SIZE(0384K) AUTOOBL( NONE ) 

•OPTIONS IN EFFECTMSOURCE EBCDIC NOLIST NOOECK NOOBJECT NOMAP NOFORMAT GOSTMT NOXREF NOALC NOANSF TERM IBM FLAG(I) 
•STATISTICS* SOURCE STATEMENTS = 76, PROGRAM SIZE = 2324, SUBPROGRAM NAME = PRINT 

""•STATISTICS* NO DIAGNOSTICS GENERATED 

end OF COMPILATION ****** 266K BYTES OF CORE NOT USED 
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LEVEL 2.3.0 (JUNE 78) 

REQUESTED OPTIONS? NOOBJ»TERMnNOXREF 


OS/36 0 FORTRAN H EXTENDED 
, ,NW1AP, » »NAME(rtAlN) .ADC NONE ) .OPTt 0) 


DATE 82.2A6/12.22.01 
. , ,FU6( I ) »SIZE( 384K ) , LC(60 ) • 


PAGE 


OPTIONS IN EFFECT 


NAME (MAIN) NOOPTiniZE 
SOURCE EBCDIC NOLIST 


LINEC0UNT(60) SIZE(038^K) AUTOOBL( NONE ) 

unnn iprT KinMA.P NOFORttAT GOSTT1T NOXREF 


NOALC NOANSF 


TERM IBM FLAG(I) 


ISN 0002 
I5N 0003 
ISN 0004 
ISN 0005 

ISN 0006 
ISN 0007 


ISN 0008 
ISN 0009 
ISN 0010 
ISN 0012 
ISN 0013 
ISN 0015 
ISN 0016 

ISN 0017 
ISN 0018 
ISN 0019 


SUBROUTINE PLOT 
IMPLICIT REAL»8(A-H,0-Z) 

OIMElwlON^ipLOTf Iz I . PLTDATI 100,201. IT< 1<W I ,RAH( <k> ,ICC 10 ) , 
1 IMAG4(5151) 

COMMON /PLOTT/ PLTDAT.IPLOT.NP 
DATA ITCl)/0/,RAH/4*0./,IC(X)/lH»/ 


C 

C 

C 


100 

no 


120 


CALCULATE NUMBER OF VARIABLES TO BE PLOTTED 


NV=0 

DO 100 1=1.42 

IF(IPLOTd) .EQ. 0) GO TO 110 


NV=NV*1 

IF(NV .EQ. 0) RETLWH 
DO 120 IP=liNV 
CALL USPLTCPLTDATtl.l 


»PLTDAT(1#IP*1 ) ,100.NP»lil»ITiRAN.ICil ► 


1 IMA64.IER) 
CONTINUE 
RETURN 
END 


00103700 

00X03800 

00103900 

00104000 

00104100 

00104200 

00104300 

00104400 

00104500 

00104600 

00104700 

00104800 

00104900 

00105000 

00105100 

00105200 

00105300 

00105400 

00105500 

00105600 

00105700 


•OPTIONS IN EFFECT*NAME(MAIN) NOOPTIMIZE 
•OPTIONS IN EFFECT-SOURCE EBCDIC NOLIST 
•STATISTICS* SOURCE STATEMENTS » 


LINECOUNT(60) SIZE(0384K> AUTODBL( NONE ) 

MOO'ECK NOOBJECT NOMAP NOFORMAT 60STMT NOXREF NOALC 
18. PROGRAM SIZE = 21736, SUBPROGRAM NAME = 


NOANSF TERM 
PLOT 


IBM FLAG(l) 


*5j^fISJICS» NO DIAGNOSTICS GENERATED 
•*••*• END OF COMPILATION <hhhhi* 
•STATISTICS* NO DIAGNOSTICS THIS STEP 


276K BYTES OF CORE NOT USED 
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APPENDIX B 


SAMPLE INPUT DATA 


This appendix provides an illustrative example of the program 
NAMELIST input data corresponding to the vehicle in Figxare 5. The 
vector geometry and inertia matrix for that vehicle are 


S 


^1 

r 


->2 

r 


-^3 

r 


■>4 

r 


tv 


b . ^ b . . ^ V 

2 ^''' 2 ^ 2 ^ 

- I i, + lb 2L)i^ + I 

b . . h , 

“2 -4 "^-^4 ‘*’ 2 ^ 

b . ^ . b , 

0 0 

(b^ t h^) 0 

0 2b^ 


12 


2 2 
(b + h ) 


0 

0 


R = 


R i^ + R 
X —1 y 


j + R k 
-^1 z —1 
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Figure 5. Example vehicle. 


The namelist input items given below correspond to the following 
dimensions and masses 

= 5 slugs 

mi = m 3 = 1 slug 
m = m. = 0.5 slug 



Also two constrained modes are to be used in the simulation. The initial 
conditions on the kinematic variables axe 


@ t 


0 


R = 1 • 10^ ft 

X 

R = 2 • 10^ ft 

Y 

R = 3 • 10^ ft 

z 

u = 0 ft/s 

V = 5 ft/s 

w = 0 ft/s 


0^ = 5 degrees 

9^ = 20 degrees 

6^ = 0 degrees 

03^ = 0 deg/s 

03^ = 10 deg/s 


0 )^ = 0 deg/s 


Note that the program in Appendix A sets the initial particle deflections, 
modal coordinates and the respective time derivatives to zero (see sub- 
routine INITL) . 

The numerical integration is to proceed from time = 0 (set in- 
ternally, see main program) to a final time of 60 seconds using an 
integration time step of 0.01 second. The print time step is to be 
6 seconds and the plot time step 0.6 second. 

The following variables are to be plotted versus time: R , 9^/ 

4 4 4 

q^, q^, V, 0)^ • 
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NAMELIST Input Data 


&INPUT M0 = 5.0, N = 4 , MASS 1.0, 0.5, 1.0, 0.5, 
RM = -0.5, 11.0, 1.0, -0.5, 21.0, 1.0, -0.5, -10. C 

- 20 . 0 , 1 . 0 , 

10 = 2.083, 0.0, 0.0, 0.0, 2.083, 0.0, 0.0, 0.0, ( 
S = -0.5, 0.5, 1.0, NT = 2 &END 

skin R = 1.E3, 2.E3, 3.E3, THETA = 5.0, 20.0, 0.0 
UVW = 0.0, 5.0, 0.0, OMEGA = 0.0, 10.0, 0.0 SEND 

S.RUN DT = 0.01, TSTOP = 60.0, DTP = 6.0, DTG = 0. 

SPLT IPLOT = 2, 5, 16, 17, 18, 20, 23 SEND 


, 1.0, -0.5, 
.833, 

> &END 
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